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1.0  INTRODUCTION 


The  series  of  programs,  known  generally  by  the  name  STAGS,  has  been  under 
development  at  the  Lockheed  Palo  Alto  Research  Laboratory  for  roughly  fifteen 
years.  The  first  version  of  STAGS  was  operational  in  about  1967  and  was  a 
finite  difference  based  program  for  the  nonlinear  analysis  of  cylindrical 
shells  with  cutouts.  This  initial  development  was  sponsored  by  LMSC 
(Lockheed).  It  was  followed  in  about  1968  by  a  special  linear  version 
restricted  to  shells  of  revolution.  Buckling  and  thermal  effects  were  added 
in  1970  followed  by  inelastic  capability,  some  finite  elements  and  more 
general  shell  geometry  by  1972.  These  programs  were  funded  by  a  number  of 
government  agencies,  but  all  went  by  the  name  STAGS.  A  new  version,  STAGS? 
(ca.  1974),  included  transient  response,  dynamic  buckling  and  could  analyze 
branched  or  segmented  shells.  The  first  version  to  be  used  by  the  structural 
analysis  community,  STAGSA,  was  released  in  1973  and  included  dynamic 
eigenvalue  analysis.  The  final  version  to  be  based  on  the  finite  difference 
method  was  released  in  1976  (STAGSC).  Since  that  time,  the  program  has 
undergone  a  major  re-writing  and  the  latest  version,  STAGSC-1,  was  released  in 
1979.  This  version  of  the  program  is  now  entirely  based  on  finite  elements 
and  future  development  will  be  along  the  same  lines.  The  development  since 
the  earlier  STAGS  version  (STAGS2  and  on)  has  been  under  the  sponsorship  of 
NASA,  Langley.  Table  1.1  summarizes  this  brief  historical  survey. 

The  present  report  deals  only  with  the  finite  element  version,  STAGSC-1.  The 
evaluation  study,  described  in  the  body  of  this  report,  follows  the  general 
methodology  described  by  Nickell  [1]  and  is  intended  to  provide  a  potential 
user  of  STAGSC-1  with  an  in-depth  description  and  critique  of  the  program.  On 
the  other  hand,  it  is  not  intended  to  be  a  "consumer  report"  by  rating  the 
program  against  other  similar  programs.  In  fact,  this  would  probably  be  a 
very  difficult,  if  not  impossible,  task  since  STAGSC-1  has  many  unique 
features  which  set  it  apart  from  most  other  finite  element  programs. 
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To  begin  with,  STAGSC- I  is  a  special  purpose  program  in  the  respect  that  its 
primary  function  is  shell  analysis  (although  a  structure  composed  entirely  of 
beams  can  also  bo  analyzed).  Secondly,  it  is  fundamentally  developed  for 
nonlinear  (geometric  and  inelastic  analysis.  Thirdly,  bifurcation  buckling  and 
dynamic  analysis  are  very  strong  features  in  STAGSC-1.  bearing  this  in  mind, 
it  is  to  be  anticipated  that  the  majority  of  users  will  be  interested  mainly 
in  the  more  sophisticated  and  less  common  analysis  problems  and  will  also  be 
aware  of  the  capabilities  of  other  candidate  programs  in  these  areas.  For 
this  reason,  the  evaluation  study  performed  for  STAGSC-1  has  done  little 
comparison  with  other  programs  and  has  concentrated  on  the  program  features, 
algorithms,  structure  and  performance. 

The  specific  tasks  which  were  performed  can  be  summarized  as  follows: 

A.  Review  of  program  documentation  (theoretical  and  users  manuals) 

G.  Program  architecture  description 

C.  Program  functional  description 

D.  Advanced  evaluation  topics; 

(1)  mesh  convergence 

(2)  eigenvalue  extraction 

(3)  transient  integration 

(4)  large  scale  nonlinear  collapse 

The  advanced  evaluation  provides  some  direct  evidence  ol  the  performance  of 
STAGSC-1  for  shell  structural  models  which  vary  in  complexity  from  the 
simplest  simile  element  model  up  to  models  with  more  I  ban  400U  degrees  of 
freedom  for  the  study  on  nonlinear  collapse.  The  mesh  convergence  study  was 
performed  with  models  up  to  630  d.o.f.,  while  the  eigenvalue  and  transient 
integration  studies  used  up  to  about  1500  d.o.f.  Not  all  capabilities  which 
are  available  were  exercised  but  the  objective  was  to  examine,  in  depth,  those 
features  which  it  is  felt  are  of  most  significance  for  the  program  in 
general.  No  specific  study  was  made  on  the  performance  of  the  modified 
Newton-Raphson  nonlinear  equation  solver.  The  reason  for  this  is  simply  that 
since  the  solver  is  used  in  most  of  the  problems  analyzed,  sufficient  evidence 
of  its  performance  would  be  generated  automatically. 
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The  criteria  for  evaluating  a  structural  analysis  computer  program  will  depend 
to  some  extent  on  the  nature  of  the  program  itself  and  hence  on  the  community 
of  potential  users.  In  the  case  of  STAGSC-1,  as  has  already  been  pointed  out, 
the  user  community  will  be  relatively  sophisticated  and  will  expect  more  than 
a  "black  box"  capability.  Thus,  program  and  theory  documentation  must  be 
accorded  substantial  weighting.  Of  course,  the  most  important  considerations 
are  still  the  ability  of  the  program  to  perform  the  types  of  analysis  for 
which  it  was  written,  accurately  and  as  economically  as  possible.  For  a 
program  which  performs  mainly  nonlinear  analysis,  relative  economy  is 
particularly  important,  since  such  analyses  tend  to  be  expensive  in  any  case. 
Ease  of  input  is  always  desirable  but  for  a  program  which  is  less  routinely 
used  it  is  not  an  overriding  factor.  More  important  for  nonlinear  analysis  is 
the  ability  to  post-process  the  results  with  as  much  freedom  as  possible. 

This  report  is  organized  into  seven  major  sections  plus  an  Appendix.  Section 
2  discusses  the  documentation  of  STAGSC-1;  Sections  3  and  4  describe  the 
program  architecture  and  functions;  Section  5  describes  program  verification 
and  Section  6  presents  the  advanced  evaluation  studies.  Conclusions  and 
recommendations  are  presented  in  Section  7.  The  Appendix  contains  reference 
diagrams  for  the  program  structure  and  also  details  of  the  element  stiffness 
matrix  modal  energy  spectrum  method  used  in  the  evaluation  of  element 
convergence. 

This  evaluation  was  performed  as  part  of  the  ISEG  program  [1]  under  ONR 
Contract  No.  N00014-79-C0825  with  Westinghouse  Advanced  Reactors  Division. 
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TABLE  1.1 
STAGS  DEVELOPMENT 


Version 
and  Date 

General  and  Added 

Capabil ities 

Sponsor 

STAGS 
(ca.  1967) 

Nonlinear  Finite  Difference  Code 
for  Cylindrical  Shells  with  Cut- 
Outs  . 

LMSC 

STAGS 

(ca.  1968/9) 

Linear  Version  for  Shells  of 
Revolution. 

NSRDC 

STAGS 
(ca.  1970) 

Bifurcation  Buckling;  Thermal 
Effects . 

SAMSO 

STAGS 

(ca.  1970/2) 

Inelastic  Analysis;  Finite 
Elements;  Extension  to  More 
General  Shell  Properties  (Grid 
Spacing) . 

AFFDL 

STAGS2 
(ca.  1972/4) 

Transient  Response;  Dynamic 
Bucklinq;  Branched  and  Segmented 
Shells. 

NASA 

Langley 

STAGSA 

1973 

Dynamic  Eigenvalue  Analysis. 

NASA 

Langley 

STAGSC 

1976 

Improved  Convergence  with 

Gridsize  for  Nonlinear 

Analysis . 

NASA 

Langley 

STAGSC 1 

1979 

Completely  Revised  Input; 

F.E.  Library  Updated  to  Include 

NASA 

Langley 

Springs,  Beams,  Shells 
(Triangle,  Quadrilateral  and 
Transition) . 


Plotting  Available  for  Geometry, 
Deformations,  Stresses,  Strains, 
etc . 


Comments 


Initial  Version. 


First  Introduction 
of  Finite  Elements. 


Versions  for  CDC 
6600  and  UNIVAC  1108 


Used  by  Structural 
Analysis  Community. 

Increased  Inaccuracy 
Due  to  Rigid  Body 
Motions. 

Program  now  Entirely 
Finite  Element 
Based. 
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2.0  PROGRAM  DOCUMENTATION 


The  documentation  pertaining  directly  to  STAGSC-1  is  at  present  confined  to  a 
users  manual  [2]  describing  the  input  and  strategy  required  to  execute  a  given 
problem.  A  theoretical  manual  also  exists  [3]  which  is  not  specific  to  the 
STAGSC-1  version.  A  third  manual,  of  example  cases  for  STAGSC,  is  available 
in  draft  form  and  is  described  as  preliminary  and  incomplete.  No  programmer's 
manual  has  so  far  been  written. 

Therefore,  the  situation  with  respect  to  documentation  appears,  at  least 
superficially,  to  be  not  entirely  satisfactory.  However,  this  is  in  part  made 
up  for  by  the  overall  high  quality  of  the  documentation  which  is  available.  A 
user,  quite  unfamiliar  with  STAGSC-1,  can  progress  to  the  point  of  successful 
execution  of  a  problem  on  the  basis  of  the  user  instructions  (Vol.  II)  alone 
[2].  Moreover,  he  can  do  this  with  some  fair  understanding  of  the  basis  of 
the  program  provided  that  he  is  reasonably  knowledgeable  on  the  subject  of 
nonlinear  finite  element  analysis.  This  is  because  there  are  sections  in  this 
volume  which  deal  with  the  important  questions  of  modeling  and  solution 
strategy.  The  theoretical  manual  was  written  when  STAGS  (STAGSC)  was  based  on 
finite  difference  theory  and  therefore  a  significant  part  of  its  content 
(20-25%)  is  no  longer  applicable.  It  clearly  needs  to  be  extensively 
rewritten  to  bring  it  up  to  date  with  the  program  but,  nevertheless,  the 
greater  part  of  it  still  provides  valuable  insight  into  the  content  and 
philosophy  of  STAGSC-1.  The  draft  of  example  problems  is  not  useful  for 
STAGSC-1  since  the  input  of  the  problems  described  was  written  for  the 
previous  version  STAGSC  and  can  no  longer  be  used.  This  section  will 
describe  and  comment  on  the  user's  and  theoretical  manuals. 

2.1  STAGSC-1  USER  INSTRUCTIONS  MANUAL 

A',  a  1  ready  noted,  the  ir.er  input  manual  is  largely  self-contained  and  provides 
the  analyst  with  explicit  instructions  for  problem  execution  together  with 
enough  background  material  to  provide  a  reasonable  understanding  of  the 
theoretical  basis  for  the  program. 
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The  manual  itself  has  ten  sections  of  which  two  minor  ones  are,  as  yet,  riot 
issued.  After  a  short  introductory  section  there  is  a  section  describing  the 
general  capabilities  of  STAGSC-1  together  witii  some  basic  concepts  on  the 
description  of  shell  surfaces.  The  third  section  provides  a  detailed 
description  of  each  input  data  card  and  groups  of  card-..  There  are  nine  such 
groups  as  t  ol lows : 

o  Summary  and  Control  Parameters 
o  Computational  Strategy  Parameters 
o  Data  Tables 
o  Geometry 
o  Discretization 
o  Boundary  Conditions 
o  Loads 

o  Output  Control 

o  t lenient  Unit 

All  input  is  free  format  and  allows  the  insertion  of  comments  which  is  a 
valuable  featjre  from  the  archival  viewpoint.  Each  input  variable  is  assigned 
a  name  which  is  usually  the  same  as  the  internal  variable  name  used  in  the 
program.  Branching  to  the  next  input  card  is  governed  by  values  of  variables 
previously  set  or  is  unconditional.  This  is  unambiguous  (at  least  as  far  as 
the  present  evaluation  is  concerned--not  all  paths  have  been  investigated)  and 
makes  for  reasonably  trouble  free  input  provided  that  the  user  takes 
sufficient  care.  At  the  end  of  each  input  card  description  there  are  a  number 
of  conditional  "go  to's"  which  lead  to  the  next  card  or  card  group.  This  kind 

of  programming  logic  for  input  preparation  is  somewhat  unusual  but  in  the 

opinion  of  the  reviewer  has  much  to  recommend  it.  At  this  point  it  is 

appropriate  to  mention  a  special  concept  in  defining  shell  geometry  which  is  a 

basic  feature  of  STAGSC-1.  A  shell  structure  may  often  be  conveniently 
described  in  terms  of  one  or  more  distinct  types  of  surface  geometry  ( e . g . , 
cylinder,  cone,  torus,  etc.).  Each  such  type  can  be  defined  in  STAGSC-1  as  a 
"shell  unit"  with  its  own  local  coordinate  system.  Shell  units  may  then  be 
connected  to  form  the  complete  structure.  Each  shell  unit  has  its  own  surface 
coordinate  grid  (rows  and  columns)  the  nodes  of  which  are  employed 
(selectively  or  otherwise)  to  define  element  connectivity.  If  a  region  of  the 
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shell  has  a  geometry  which  does  not  lend  itself  to  this  treatment,  nodes  and 
elements  may  be  specified  directly  to  form,  what  is  referred  to  in  STA6SC-1 
terminology,  as  an  element  unit.  Further  details  are  discussed  in  Section  4.2. 


The  summary  and  control  parameter  group  provide  a  title  and  define  the 
analysis  type.  The  next  group,  computational  strategy  parameters,  define  the 
following:  load  incrementation;  eigenvalue  extraction  strategy  (spectral 
shift,  eigenvalue  range);  time  history  forcing  function  input  and  integration 
method;  basic  mesh  summary  for  each  shell  unit  (rows  and  columns);  shell  unit 
interconnections;  element  unit  summary  (element  types  and  number).  The  third 
group  (data  tables)  specifies  all  material  properties,  cross-sections  of  beam 
elements,  shell  wall  construction  (multilayer,  composite,  etc.)  and  a  table  of 
real  and  integer  constants  for  use  in  user  coded  subroutines  if  required. 

Shell  unit  surface  geometry  is  defined  in  the  fourth  group.  This  may  be 
chosen  from  a  library  of  standard  forms  (11  surfaces)  or  generated  by  a  user 
subroutine.  A  thirteenth  option,  to  fit  a  surface  to  defined  points  by  spline 
functions,  is  indicated  but  not  currently  available.  Shell  wall  type  is  also 
specified  in  this  group  together  with  the  type  of  strain-displacement 
relationships  (linear  or  nonlinear),  elastic  or  elastic-plastic  material 
behavior  and  initial  imperfections  (if  of  trigonometric  form). 

The  next  major  input  group  controls  the  finite  element  discretization.  The 
analyst  may  use  the  mesh  defined  in  the  shell  geometry  or  specify  an  overlying 
mesh  of  elements  which  picks  shell  mesh  nodes  selectively  for  element 
connectivity.  Patches  of  elements  can  be  defined  which  allow  changes  in 
element  type;  also  segments  of  the  mesh  can  be  defined  with  different 
spacing.  In  addition,  the  user  may  define  the  element  connectivity  through  a 
user  subroutine.  Also  specified  in  this  group  are  the  location  of  discrete 
stiffeners  which  may  be  eccentric  with  respect  to  the  shell  surface  and  can  be 
skewed  with  respect  to  the  shell  surface  mesh. 

Boundary  conditions  for  shell  units  are  specified  in  the  sixth  group.  These 
can  be  of  a  standard  type  (simple  support,  clamped,  etc.)  along  a  specified 
edge  of  the  shell  or  they  may  be  selected  degree  of  freedom  types  along  the 
edge  which  can  be  designated  fixed  or  free.  An  additional  distinction  may  be 
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drawn  between  boundary  conditions  applicable  to  incremental  displacements  or 
basic  displacements  (pre-buckl ing,  prestress).  The  seventh  group  controls  the 
applied  loads,  displacements  and  initial  conditions  (for  dynamic  analysis).  A 
concept  is  introduced  here  which  is  unusual  in  structural  analysis  programs, 
i.e.,  the  specification  of  loads,  displacements  or  initial  conditions  in  two 
categories,  A  and  B,  which  are  independent.  They  are  scaled  by  two 
independent  load  factors  and  Pg.  The  purpose  of  this  is  specifically 
the  determination  of  bifurcation  buckling  behavior  for  systems  where  there  is 
a  fixed  load,  or  prestress,  defined  by  System  B  and  an  increasing  load  defined 
by  System  A.  At  bifurcation  the  total  stress  is  obtained  from  the  sum  of 
System  B  stresses  and  System  A  multiplied  by  the  eigenvalue.  For  non-buckling 
analysis  it  is  not  necessary  to  specify  both  A  and  B  systems.  For  a  new  user, 
reading  the  documentation  unaided,  this  concept  can  be  rather  puzzling  at 
first  and  the  documentation  should  provide  a  little  more  discussion  than  it 
does.  However,  in  all  other  respects  the  description  of  loads  input  is  guite 
satisfactory. 

Output  control  is  input  on  the  eighth  group  of  data  cards.  Basically,  the  user 
can  control  the  frequency  of  printing  displacements,  stress  resultants, 
strains,  stresses  and  point  forces.  The  frequencies  for  each  of  the  above 
quantities  are  independently  specified.  These  controls  govern  printout  for 
all  elements  in  a  given  unit  and  must  be  specified  separately  for  each  unit. 

In  addition,  certain  selected  stress  or  displacement  components  may  be  printed 
at  every  step.  Apart  from  this  control,  the  manual  recommends  the  use  of  the 
post-processor  STAPL.  This,  however,  is  geared  mainly  to  the  generation  of 
contour  plots  and  does  not  currently  fulfill  many  post-processing  needs,  e.g., 
time  history  of  stresses  or  displacements  or  spatial  distubution  of  stresses 
or  displacements.  A  significant  aid  in  this  respect  would  be  detailed 
descriptions  of  the  model  and  solution  data  files  (M0D  and  S0D)  so  that  the 
data  may  be  processed  according  to  the  user's  needs  in  a  separate  program. 

Such  descriptions  are  not  given  in  the  documentation.  All  of  the  preceeding 
input  pertaining  to  the  shell  unit  must  be  repeated  for  each  shell  unit  in  the 
structure.  Each  shell  unit,  being  independently  defined,  may  be  totally 
different  with  regard  to  each  section  of  the  input  data  which  is  part  of  the 
shell  unit  (geometry,  elements,  wall,  material,  etc.) 
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The  final  input,  group  is  for  definition  of  an  element  unit.  This  group 
specifies  user  point  coordinates,  element  connectivity,  point  forces  at  nodes 
and  distributed  forces  on  elements.  This  part  of  the  input  is  rather 
unsophisticated  since  input  by  cards  requires  specification  of  each  and  every 
node  point  and  its  coordinates  arid  each  element  and  its  nodes.  Ihr  only 
alternative  to  this  is  to  define  node  geometry  and  element  connectivity  by 
means  of  user  subroutines. 

The  fourth  section  of  the  manual  describes  the  function  and  format  of  each  of 
the  thirteen  available  user  subroutines.  The  descriptions  are  reasonably 
complete  and  each  one  provides  a  coded  example. 

Section  five  describes  the  input  for  the  post-processor  S1APL.  lhe  input 
required  is  in  free  format,  as  in  STAGS,  and  is  straightforward.  Currently 
irioperational  features  are  marked  by  an  asterisk. 

Section  6,  dealing  with  modeling  and  strategy,  is  one  of  the  most  important 
sections  in  the  manual.  It  provides  the  user  with  valuable  background 
information  with  regard  to  the  type  of  analysis  which  should  be  performed  to 
solve  a  particular  problem.  It  is  in  this  section  that  one  may  discern  the 
basic  philosophy  of  the  STAGS  series  of  programs.  Put  concisely,  the 
fundamental  viewpoint  is  that  the  purpose  of  shell  analysis  is  determination 
of  the  structural  failure  modes.  For  thin  shells,  this  usually  means  either 
static  or  dynamic  collapse  of  the  shell  wall  rather  than  a  straightforward 
exceeding  of  material  stress  or  strain  limits.  Thus,  S1AGSC-1  is  the  result 
of  an  evolution  of  a  basic  nonlinear  shell  analysis  capability  with  strong 
emphasis  on  bifurcation  buckling.  The  dynamic  counterpat  L  has  also  been 
developed  and  plastic  material  behavior  included. 

An  indication  of  the  degree  of  development  of  the  nonlinear  capability  is 
provided  by  the  statement  (in  Section  6)  ".  .  .  it  is  suggested  that  nonlinear 
static  analysis  be  used  as  a  matter  of  course.  .  ."  on  the  grounds  that  if 
nonlinear  effects  are  insignificant,  the  nonlinear  solution  will  converge  in 
one  or  two  iterations  and  the  extra  cost  will  be  negligible.  General  reasons 
arc  discussed  for  the  non-convergence  of  a  problem.  A  few  paragraphs  are  also 
devoted  to  the  problems  of  bifurcation  buckling  analysis. 
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The  next  subsection  discusses  mode  1  i riy  witti  STAGSC-)  using  tin-  shell  unit 
concept  for  geometry  description  and  also  the  use  of  the  capability  to 
incorporate  either  discrete  or  "smeared"  stiffeners.  Discretization  is 
discussed  next  and  is  basically  a  description  of  the  element  library  available 
in  STAGSC-1.  This  is  of  vital  importance  since  it  is  t fie  only  place  in  the 
documentation  where  the  elements  are  described.  It  is  adequate  for  general 
information  but  does  not  provide  the  depth  of  detail  which  should  be  contained 
in  <i  propel  theoretic. 1 1  desc  r  ipt  ion . 

The  final  subsection  deals  with  computational  strategy.  First,  eigenvalue 
analysis  for  both  bifurcation  buckling  and  vibration  frequencies  is 
discussed.  The  subtleties  of  bifurcation  buckling  are  highlighted  (existence 
of  negative  as  well  as  positive  eigenvalues)  and  also  the  choice  of  spectral 
shift  values.  The  general  strategy  of  nonlinear  analysis  is  covered  next. 

The  user  specified  parameters  NCUT  (total  number  of  times  the  load  increment 
may  be  halved)  and  NtWT  (total  number  of  times  the  factored  matrix  may  be 
computed)  are  also  dealt  with  in  some  detail  since  they  are  quite  subtle  in 
their  effect  on  the  modified  Newton-Raphson  procedure.  Also,  the  internal 
logic  of  the  solution  algorithm  is  quite  complex  in  the  way  it  makes  decisions 
with  regard  to  refactoring  or  cutting  the  load  step.  As  a  result,  even  given 
a  careful  reading  of  this  section,  the  user  may  feel  that  he  can  exert  more 
control  over  the  solution  procedure  than  is  actually  the  case.  The  parameters 
UELX  (tolerance  criterion  for  displacement  increments)  and  WUND  (relaxation 
factor)  are  also  user  specified  and  their  use  discussed. 

Finally,  the  use  of  the  transient  integration  operators  is  covered.  Most,  of 
the  discussion  deals  with  the  central  difference  (explicit)  operator  and  the 
question  of  its  conditional  stability  and  how  to  estimate  the  critical  time 
step.  The  discussion  of  the  implicit  operators  is  less  thorough  and  more 
guidance  on  the  appropriate  choice  for  a  given  problem  would  be  useful. 
Surprisingly,  there  is  hardly  any  reference  to  plasticity  calculations  in  this 
section.  Given  that  plasticity  is  certain  to  enter  into  most  nonlinear 
calculations  the  omission  is  serious.  It  is,  however,  true  that  there  is  not 
much  opportunity  for  the  user  to  exert  control  over  the  plasticity 
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computations  (other  than  via  the  specification  of  the  uniaxial  stress-strain 
curve)  and  that  the  program  handles  them  automatically.  Some  discussion  on 
limitations  should  be  provided. 

Section  7  of  the  manual  attempts  to  provide  user  guidance  on  the 
interpretation  of  output  from  STA6SC-1.  Special  problems  which  may  arise 
durinq  linear,  eigenvalue,  nonlinear  and  transient  analyses  are  presented  and 
appropriate  actions  suggested.  While  such  a  list  can  never  be  complete,  the 
situations  presented  are  basic  and  provide  general  guidance  on  how  the  user 
should  handle  difficulties  in  execution. 

Section  8  (Index  to  Volumes  1  and  2)  and  9  (Execution  Control)  are  in 
preparation.  An  index  to  the  users'  manual  (Vol.  2)  is  definitely  a 
necessity.  A  disadvantage  of  the  way  the  input  is  structured  is  that  it  is 
not  easy  to  locate  the  place  in  the  manual  where  a  specific  item  is 
discussed.  For  example,  in  order  to  determine  what  actions  are  required  to 
create  a  restart  file  it  is  necessary  to  trace  through  some  of  the  input 
groups  in  detail.  A  well  constructed  index  would  help  the  user  of  only 
moderate  experience  very  considerably.  The  proposed  section  on  execution 
control  is  perhaps  of  lesser  importance  since  this  must  be  different  for  each 
system  on  which  STAGSC-1  is  installed. 

Finally,  Section  10  (Minimanual)  provides  a  useful  summary  of  all  the  input 
records  and  the  associated  variable  names.  This  can  be  also  used  as  a 
stop-gap  until  a  proper  index  is  available. 

2.2  THEORETICAL  MANUAL 

As  indicated  previously,  the  currently  available  theoretical  manual  was  not 
written  in  the  context  of  the  STAGSC-1  program.  The  latest  version  of  the 
program  to  which  it  truly  applies  is  STAGSC,  which  is  a  finite  difference 
program  whereas  STAGSC-1  is  wholly  finite  element.  Thus,  there  are  parts  of 
this  volume  which  are  still  valid  and  others  which  are  irrelevant  and  it  is 
part  of  the  purpose  of  this  subsection  to  identify  the  portions  which  remain 
valid.  The  manual  is  organized  into  ten  sections  plus  an  appendix  as  follows 
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o  Introduction 

o  Summary  of  Theory 

o  Constitutive  Relations  (Elastic) 
o  The  Theory  of  Plasticity 

o  Geometric  Nonlinearities 

o  Discretization  Procedures 

o  Beams  and  Stiffeners 

o  Constraints  and  Transformations 

o  Solution  Procedures 

o  Program  Organization 

Appendix 

The  current  version  of  the  manual  is  incomplete;  Sections  7,  8,  10  and  two 
subsections  of  the  Appendix  are  not  yet  available. 

The  introductory  section  starts  off  with  some  general  observations  on 
computer-based  structural  analysis  and  goes  on  to  list  the  basic 
approximations  inherent  in  the  STAGS  program.  It  must  again  be  emphasized 
that  many  comments,  observations  and  whole  sections  of  this  volume  are  not 
appropriate  to  the  STAGSC-1  version.  This  is  already  apparent  in  the 
limitations  described  for  first  order  shell  theory.  Since  STAGSC-1  is  finite 
element  based  and  there  are  no  curved  shell  elements  available  in  the  program, 
shell  theory  cannot  be  a  part  of  the  element  formulation.* 

Section  2  provides  a  summary  of  the  basic  theoretical  principles  embodied  in 
STAGS.  First,  there  is  a  short  discussion  on  variational  principles  in  terms 
of  Hamilton's  principle  for  dynamics  which  is  shown  to  degenerate  to  the 
principle  of  minimum  potential  energy  for  static  systems.  A  point  of  interest 
is  that  the  conditions  under  which  "live"  pressure  loads  are  conservative  are 
delineated.  Thus,  it  is  possible  to  include  these  in  an  analysis  based  in  the 
variational  approach  which  is  necessarily  confined  to  conservative  systems. 

The  basic  theory  of  elasticity  constitutive  equations  are  presented  and  the 


*The  elements  are  all  flat  plate  triangles  or  quadrilateral  in  which  there  is 
no  coupling  of  membrane  and  bending  behavior  within  the  element. 


091 7B-84B : 2 
(S3034)  12 


12 


point  ni.it!*'  that  the  STAGS  formulation  is  bast'd  on  engineer  inq  strain. 

Kinematic  strain  displacement  relationships  are  developed  for  small  strains 
but  with  moderate  rotations  (<0.3  radians).  These  are  Green's  strain  in 
Laqrangian  coordinates. 

The  third  section  presents  a  development  of  generalized  elastic  constitutive 
relationships  for  combinations  of  shell  wall  and  stiffeners. 

Section  4  provides  a  reasonably  up-to-date  discussion  of  plasticity.  First, 
the  classical  development  of  plasticity  theory  is  considered,  incremental  and 
deformation  theories  are  discussed  with  deformation  theory  beinq  discarded  for 
general  non-monotonic  loading.  Reverse  plasticity  and  the  Bauschinger  effect 
are  also  treated.  On  the  basis  of  improved  correlations  with  experiment 
(compared  with  isotropic  or  kinematic  hardening),  the  White-Bessel ing 
(mechanical  sublayer)  theory  is  selected  as  being  representative  of  more 
modern  theories  of  plastic  work  hardening  behavior.  It  is,  however,  pointed 
out  that  for  complex  loading  (non-monotonic  and  non-proportional),  little  is 
known  about  the  applicability  of  currently  available  theories.  The  solution 
of  plasticity  problems  using  the  pseudo-force  and  the  tangential  stiffness 
methods  is  discussed.  It  is  also  indicated  that  STAGS  uses  a  combined 
approach  by  normally  using  the  pseudo-force  technique  and  updating  the 
stiffness  matrix  for  plasticity  effects  whenever  convergence  behavior 
indicates  reformulation  of  the  stiffness  is  required.  This,  however,  is  quite 
misleading  since  the  STAGSC-1  program  does  not  include  plasticity  corrections 
when  formulating  the  stiffness.  An  additional  feature  of  the  plasticity 
computations  is  the  so-called  subincrement  approach.  This  is  discussed  later 
in  this  report  in  Section  4.5  of  the  functional  description. 

Section  5  is  entitled  "Geometric  Nonlinearities"  but  concentrates  entirely  on 
the  question  of  structural  stability.  Nevertheless,  the  discussion  is  very 
thorough  and  reveals  where  the  greatest  depth  of  expertise,  applied  to  the 
development  of  STAGS,  lies.  The  concept  of  structural  stability  is  developed 
in  terms  of  primary  and  secondary  loading  paths  and  the  bifurcation  point. 
Various  criteria  for  instability  are  discussed  followed  by  the  consequences  of 
instability  and  pre-  and  post-buckling  behavior. 
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The  rest  of  the  section  discusses  static  stability  analysis  with  a  large 
amount  of  qualitative  detail.  A  short  historical  survey  of  analysis  methods 
is  followed  by  a  dissertation  on  flat  plate  buckling  including  post-buckling 
behavior.  This  is  followed  by  analysis  of  shells  of  revolution  and  general 
shell  buckling.  For  the  user  who  is  mainly  preoccupied  with  shell  stability 
problems  this  section  provides  an  excellent  background. 

The  next  section  (Section  6)  is  a  major  discussion  on  discretization 
procedures.  Unfortunately,  much  of  it  is  no  longer  of  direct  consequence  to 
STAGSC-1  since  it  discusses  finite  difference  procedures  in  some  depth. 

Also,  the  discussion  on  finite  elements  is  somewhat  out-of-date  with  respect 
to  the  specific  elements  available  in  STAGSC-1.  Nevertheless,  the  material 
presented  is  of  high  quality  and  the  contents  of  the  section  will  be  reviewed 
here,  at  least  in  part. 

The  first  three  subsections  give  an  overview  of  standard  methods  for  numerical 
differentiation  and  integration.  This  is  followed  by  a  discussion  in  some 
depth  of  numerical  solution  procedures  with  considerable  emphasis  on  finite 
difference  methods.  Finite  element  procedures  in  general  are  discussed  mainly 
from  the  standpoint  of  continuity  requirements  for  convergence.  Topics 
covered  are  CQ,  continuity  requirements,  conforming  and  non-conforming 
elements,  order  of  convergence  and  the  patch  test.  The  presentation  is 
general  and  not  specific  to  the  STAGSC-1  program.  With  regard  to  the  finite 
difference  versus  the  finite  element  approach,  the  point  of  view  of  the 
developers  is,  to  quote  verbatim; 

"There  is  no  clear  distinction  between  the  finite  element  method  and  the 
finite  difference  energy  method.  It  seems  reasonable  to  define  as  a 
finite  element  method  a  discretization  scheme  in  which  the  displacement 
pattern  inside  the  element  is  determined  without  the  use  of  nodal 
freedoms  outside  the  closed  domain  of  the  element." 

In  the  opinion  of  the  reviewer,  this  perspective  sheds  some  light  on  the 
reasons  for  the  choice  of  program  architecture  in  STAGSC-1.  This  will  be 
described  in  some  detail  in  Section  3  but  it  can  be  said  at  this  point  that  it 
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is  quite  different  from  conventional  finite  element  programs.  Under  the 
heading  of  "Special  Problems",  the  following  topics  are  discussed  briefly; 
effects  of  reduced  integration  in  producing  spurious  mechanisms;  strains 
produced  by  rigid  body  motion  ( super-parametric  elements);  the  undesirability 
of  convergence  from  "below"  for  buckling  analysis. 

Finite  difference  schemes  and  some  finite  elements  are  discussed  next.  A 
so-called  "curved"  finite  element  (STA6C)  with  incompatible  ("bubble")  modes 
is  discussed  but  does  not  appear  to  have  been  included  in  any  available 
version  of  the  program.  A  flat  plate  quadrilateral  element  (STAGF)  is 
introduced  next  and  the  compatibility  problems  associated  with  modeling  curved 
shell  surfaces  with  flat  elements  is  discussed  in  some  detail.  It  is 
concluded  that  for  displacement  continuity,  cubic  variation  of  in-plane 
displacement  components  normal  to  an  edge  is  required  in  order  to  match  the 
cubic  variation  of  transverse  (bending)  displacements.  In  addition,  degrees 
of  freedom  corresponding  to  average  in-plane  shear  strain  and  rotation  about 
the  surface  normal  are  required  at  corner  nodes.  Thus,  the  complete  element 
is  specified  with  seven  freedoms  at  corner  nodes  and  four  at  midside  nodes  (32 
total).  No  derivation  is  provided. 

The  merits  of  Ahmad  type  elements  for  thin  shell  analysis  are  discussed  and 
finally  the  Clough-Fel ippa  triangular  and  quadrilateral  elements.  The  section 
concludes  by  listing  elements  to  be  included  in  STAGS  as  follows: 

o  Flat  quadrilateral  STAGF 

o  Ahmad  type  elements 

o  Clough-Fel ippa  triangle  and  quadrilateral 

As  will  be  detailed  in  Section  4,  these  are  not  the  elements  that  are  in  the 
STAGSC-1  program,  but  it  is  clear  that  they  are  derivative  versions  of  the 
Clough-Fel ippa  series.  Ahmad  elements  are  not  yet  available. 

Section  7  (Reams  and  Stiffeners)  and  8  (Constraints  and  Transformations)  are 
described  as  being  in  preparation. 
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bee  t  ion  9,  entitled  "Solution  Procedures",  is  again  one  of  L  fie  major  parts  of 
the  theoretical  manual.  Topics  discussed  are  as  follows: 


o  Expressions  for  strain,  kinetic  and  total  potential  energies 
o  Linear  equation  systems 
o  Nonlinear  equation  systems 
o  Eigenvalue  analysis 
o  Transient  analysis 

This  section  is  probably  the  most  mathemat ical ly  detailed  part  of-the 
theoretical  manual;  it  begins  with  statements  of  strain,  kinetic  and  potential 
energies  and  introduces  the  concept  of  a  stiffness  operator  L(x)  which  is 
defined  ,is  the  first,  variation  of  the  total  potential  energy.  L ( x )  is  defined 
as  generally  nonlinear,  thus,  there  is  a  departure  from  the  more  familiar 
ideas  of  linear  and  nonlinear  stiffness  matrices  in  the  subsequent 
developments.  The  solution  of  linear  systems  of  equations  is  presented  in 
terms  of  conventional  triangular  decomposition  followed  by  forward  arid 
backward  substitution.  There  is  also  an  important  discussion  of  problems 
encountered  in  the  solution  process  due  to  ill-conditioning.  This  gives 
valuable  insight  into  the  meaning  of  the  various  diagnostic  messages  which  may 
be  output  during  STAGS  execution. 

I  tie  solution  of  nonlinear  equation  systems  is  also  discussed  in  depth, 
beginning  with  a  discussion  of  the  relative  merits  of  regular  and  modified 
Newtori-kaphson ;  successive  substitution  with  nonl  inear ities  on  t fie  right  hand 
side  only;  tangent  stiffness  incremental  method  with  residual  load  correction; 
dynamic  relaxation  and,  finally,  energy  search  methods.  It  is  concluded  that 
the  Newton-kaphson  methods  include  the  tangent  stiffness  methods  as  special 
cases  and  that  dynamic  relaxation  is  not  competitive.  The  discussion  then 
goes  on  to  develop  equations  for  regular  and  modified  Newton-kaphson  in  terms 
of  the  operator  L(x)  (there  are,  unfortunately,  two  errors  in  the  equations 
which  require  correction).  It  is  stated  that  the  user  can  involve  either 
regular  or  modified  Newton-Raphson  but  this  does  not  appear  to  be  operational 
in  the  version  evaluated  (see  Section  4.5).  It  is  also  stated  that  the  user 
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can  choose  to  include  material  nonlinearities  in  reformulation  of  the 
stiffness  but,  again,  this  is  not  available  in  STAGSC-1.  Plastic 
nonlinearities  are  included  as  pseudo-force  contributions  to  the  loads  vector. 

Eigenvalue  analysis  is  introduced  in  very  general  terms  through  use  of  the 
nonlinear  stiffness  operator  L(x).  While  concise,  this  treatment  tends  to  be 
somewhat  obscure  to  the  less  mathematically  inclined  structural  analyst. 

Having  stated  the  linear  eigenvalue  problem,  the  generation  of  the  associated 
matrices  is  developed  in  terms  of  second  derivatives  ("second  variation")  of 
the  potential  energy.  This  step  is  revealing  with  respect  to  the  programming 
of  STAGS  since  it  appears  that  this  is  the  basis  of  the  algorithms 
implemented.  The  solution  of  the  eigenvalue  problem  is  described  in  terms  of 
the  inverse  power  method  including  a  spectral  shift.  However,  the  actual 
method  employed  in  STAGS  is  not  described  in  detail.  This  is  unfortunate, 
since  the  detailed  treatment,  which  was  to  be  included  in  the  Appendix  to  the 
theory  manual,  has  not  been  written.  As  is  described  in  Section  4.7,  the 
actual  method  used  is  a  variation  on  the  subspace  iteration  method,  which  is  a 
relatively  recent  development. 

The  final  topic  is  transient  integration,  which  is  one  of  the  strong  features 
of  STAGSC-1.  The  treatment  is  good  and  is  substantially  more  clear  than  some 
of  the  preceding  sections.  Again,  the  subject  is  introduced  in  general  terms 
with  a  discussion  of  the  solution  of  initial  and  boundary  value  problems. 
Explicit  and  implicit  integration  methods  are  defined.  The  explicit  (central 
difference)  algorithm  is  presented  and  its  conditional  stability  discussed. 
There  is  also  a  lengthy  discussion  of  implicit  methods  and  their  stability. 

Data  are  presented  for  the  stability  boundaries  (applicable  to  a  linear 
analysis)  for  a  number  of  schemes,  viz.  Park,  Wilson-9,  Houbolt,  Gear's  2nd 
and  3rd  order  and  the  trapezoidal  method.  With  regard  to  their  stability  for 
nonlinear  analyses  it  is  pointed  out  that  the  criteria  are  not  exactly  valid 
and  that  no  method  exists  which  is  unconditionally  stable.*  It  may  be 
difficult,  therefore,  to  distinguish  between  physical  and  numerical 
instability  in  the  nonlinear  case.  It  is  concluded  that  energy  balance  checks 
are  advisable;  however,  STAGSC-1  is  not,  as  yet,  provided  with  this  capability. 


*See  further  discussion  in  Section  4.8. 
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The  Appendix  section  is  incomplete,  since  it  contains  only  a  section  on  shell 
theory,  which  is  not  required  for  S1AGSC-1.  There  should  also  he  more 
detailed  information  on  the  handling  of  plasticity,  element  formulations,  the 
eigerisolver  and  transient  integration  operators. 

2.3  CONCLUSIONS 

It  may  be  concluded  from  the  foregoing  that  there  are  some  serious 
deficiencies  in  the  documentation  for  STAGSC-1.  These  are: 

o  lack  of  any  program  description 
o  partial  obsolescence  of  the  theoretical  manual 
o  lack  of  a  problem  demonstration  manual 

On  the  positive  side,  it  may  be  fairly  stated  that  the  quality  of  the  existing 
manuals  is  high  and  that  they  are  written  with  the  sophisticated  user  in 
mind.  Since  STAGSC-1  is  basically  a  nonlinear  program,  this  is  the  right 
approach.  However,  the  availability  of  a  complete  set  of  manuals  cannot  be 
predicted  at  the  time  of  writing  and  it  is  to  be  hoped  that  some  greater 
priority  will  be  given  by  the  developers  to  this  highly  important  aspect  of 
the  program. 
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3.0  PROGRAM  ARCHITECTORE 


STaGSC-1  is  the  latest  version  of  the  STAGS  series  of  programs  developed  over 
the  last  13  years  or  so.  It  is,  however,  a  major  departure  from  previous 
versions  since  it  is  now  entirely  finite  element,  rather  than  finite 
difference  based.  The  program  itself  has  been  completely  rewritten  and 
probably  bears  little  resemblance  to  the  earlier  versions  (although  a 
comparison  of  the  current  and  earlier  source  listings  has  not  been  made  to 
verify  this).  However,  inspection  of  the  job  control  deck  required  for 
execution  shows  that  there  are  major  differences  in  organization  of  the 
program.  The  most  obvious  feature  is  that  STAGSC-1  consists  of  two  programs, 
referred  to  as  STAGS1  and  STAGS2.  STAGS1  is  a  pre-processor  which  reads  the 
input  data,  generates  the  finite  element  model,  derives  the  nodal  forces  and 
creates  a  file  to  preserve  the  data  base  for  execution.  STAGS?  then  takes 
over  arid  performs  the  execution.  A  new  feature  is  a  separate  post-processor, 
SImPL,  winch  provides  geometry  plots  of  the  model  and  contour  plots  of  the 
solution  variables.  STAPL  is  not  currently  fully  operational  and  has  been 
excluded  from  the  evaluation.  This  section,  therefore,  will  describe  in  some 
depth  the  major  features  of  the  architecture  of  STAGSC-1. 

The  performance  of  this  part  of  the  evaluation  has  been  hampered  by  the 
complete  lack  of  documentation  on  the  programming  aspects  of  STAGSC-1.  In 
addition,  the  lack  of  a  revised  theoretical  manual  has  also  made  it  difficult 
for  the  most  part  to  establish  the  algorithms  embedded  in  the  program.  This 
is  especially  true  with  respect  to  the  creation  of  the  stiffness  and  the 
solution  of  equations  since  the  procedures  involved  are  unusual. 

3.1  GENERAL  DESCRIPTION 

As  already  mentioned,  the  STAGSC-1  program  is  a  system  of  three  separate 
programs,  STAGS1,  STAGS2  and  STAPL.  STA6S1  and  STAGS2  are  normally  executed 
in  tandem  with  STAPL  following  or  run  in  a  stand-alone  mode.  The  version  of 
MAGSC-l  which  was  evaluated  was  configured  for  the  COL- /600 .  Both  STAGS1  and 
S1AGS2  are  highly  modularized  with  primary  and  secondary  overlays.  STAGS1  has 
eight  primary  and  three  secondary  overlays  while  STAGS2  has  eight  primary  and 
four  secondary. 
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The  organization  of  the  program  for  the  CDC-7600  using  the  SCOPE  2  system 
(Version  2.1.5)  utilizes  the  library  file  concept.  Three  files  are  created; 
one  for  routines  used  in  STAGS1,  one  for  STAGS2  and  the  third  for  routines 
used  in  both  STAGS1  and  STAGS2. 

The  source  coding  is  FORTRAN  with  a  few  routines  written  in  assembly  language 
(COMPASS).  However,  equivalent  FORTRAN  coding  is  supplied  in  the  form  of 
comment  statements.  The  degree  of  commenting  in  the  program  is  somewhat 
variable.  The  majority  of  routines  have  at  least  some  description  of  their 
function;  a  number  of  others  are  rather  fully  commented  while  a  few  have  no 
comments  at  all.  Overall,  the  commenting  is  adequate,  but  is  not  sufficient 
to  obviate  the  need  for  proper  programming  documentation. 

3.2  STAGS1  -  PRE-PROCESSOR 

The  overlay  structure  of  STAGS1  is  shown  schematically  in  Figures  3.1  and 
3.2.  The  (0,0)  overlay  is  a  short  main  program  which  serves  to  load  primary 
and  secondary  overlays  as  required  by  the  input  data.  As  indicated  by  Figure 
3.1,  the  main  program  calls  overlays  (1,0),  (2,0)  and  (7,0)  directly  and  also 
the  other  primary  overlays  through  the  subroutine  PREVU.  The  direct  calls  are 
for  the  purpose  of  determining  core  storage  requirements.  In  addition,  the 
main  program  saves  the  model  data  generated  on  the  file  MOT. 

Of  the  eight  primary  overlays,  three  are  always  loaded  (6,0),  (7,0)  and  (8,0); 
these  summarize  the  model  input  data,  prepare  the  element  stiffness  data  and 
save  the  data  base.  The  remaining  five  primary  overlays  are  loaded 
selectively  depending  on  whether  the  model  consists  of  shell  units,  an  element 
unit  or  both. 

Similarly,  the  secondary  overlays  are  loaded  selectively  depending  on  the 
class  of  elements  called  for  at  input  time  (beam,  quadrilateral  or  triangle). 

A  detailed  set  of  diagrams  showing  the  calls  made  to  all  subroutines  are 
provided  in  Section  9  (Appendix).  The  main  program  (Figure  9.1)  also  loads 
the  user  subroutines  (LUSER1)  which  may  be  required  for  definition  of  the 
structural  model  (Figure  9.2).  Major  subroutine  functions  are: 


09748-86B : 2 
(S3034)  2 


20 


o  PHASE! -controls  all  pre-processing 

o  LOVO-loads  overlays  from  level  0 

o  AOATA-sets  up  data  statements 

o  BDATA-sets  up  node  number  and  logical  freedom  lists 
3.2.1  OVERLAY  FUNCTIONAL  DESCRIPTIONS 

The  primary  and  secondary  overlays  will  each  be  described  briefly  in  this 
section. 

A.  Overlay  (1,0)  -  Program  OVSU  -  Shell  Unit  Generation.  Figure  9.4 
shows  the  schematic  representation  of  the  routine  calls  made  in  this 
overlay.  The  function  of  this  overlay  is  to  generate  the  underlying 
mesh  for  the  shell  units  and  the  stream  of  associated  elements. 
Control  over  the  generation  is  via  subroutine  GENSU  which  calls  the 
major  routines  GINPT,  MSHGEN,  SUN  and  SUE.  GINPT  selects  the  type  of 
shell  unit  specified  (cylinder,  torus,  etc.  or  defined  by  user 
subroutine  LAME)  and  defines  its  global  orientation,  wall 
construction  and  reference  surface  imperfections  (WIMP).  MSHGEN 
generates  the  underlying  mesh  when  the  spacing  of  the  gridlines  is 
non-uniform.  The  actual  nodal  coordinates  are  computed,  saved  in  a 
node  "table"  and  also  printed  out  by  the  subroutine  SUN;  SUN  also 
checks  for  consistency  in  the  mesh  specified  by  the  user.  Element 
connectivity  is  generated,  stored  and  checked  for  consistency  by 
SUE.  Schematics  of  calls  made  by  SUN  and  SUE  are  given  in  Figures  9.5 
and  9.6.  The  following  routines  establish  the  element  data  in  the 
shell  unit  configuration  tables: 

ATRIA  -  alternating  triangularization  for  quadrilateral  elements 

BEAM  -  beam  elements 

QUAD  -  quadri lateral  elements 

QUIN  -  transition  elements 

STIFIN  -  discrete  stiffeners 

TRI  -  triangular  elements 
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6.  Overlay  (2,0)  -  Program  OVEU  -  Element  Unit  Generation. 


The  logic  for  element  unit  generation  is  similar  to  that  for  the 
shell  unit  but  simpler  because  there  are  no  built  in  geometries 
(Figure  9.7).  Thus,  subroutine  EUN  reads  in  node  point  numbers  and 
coordinates  directly  or  through  the  user  subroutine  USRP1 . 

Subroutine  EUE  generates  element  connectivity  data  according  to 
element  type. 

C.  Overlay  (6,0)  -  Program  OVIS  -  Model  Input  Summary 

This  overlay  creates  files  which  contain  all  the  additional  data 
required  for  the  execution  of  the  analysis.  Input  data,  which  is 
entered  in  free-format  is  interpreted  by  subroutine  CARDS  ana 
translated  into  an  internal  format.  The  master  subroutine,  PREMIS, 
assembles  the  input  data  defining  the  analysis  type,  loading  and 
solution  strategies,  structural  model  definition  and  provides 
descriptive  output.  In  addition,  PREMIS  creates  a  beam  cross-section 
properties  file,  material  properties  file  and  shell-wall  construction 
file.  Data  defining  the  method  of  time  integration  (if  used)  is  also 
loaded  and  saved  (LOADT).  The  following  subroutines  perform  major 
functions  as  follows: 

CARDS  -  interprets  free-form  input 

ESP1D  -  generates  White-Bessel ing  plasticity  data 

RCONST  -  reads  constraint  conditions 

TAB  -  tabulates  beam  section  properties 

TAM  -  tabulates  user  materials 

TAP  -  tabulates  user  parameters 

TAW  -  tabulates  shell  wall  properties 

Figure  9.8  shows  the  subroutine  calls  for  this  overlay. 
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I).  Overlay  (7,0)  -  Proqram  OVVU  -  Element  Stiffness  Preprocessor 

The  descriptive  comments  for  this  overlay  use  the  terminology  “unit 
prevar iational  overlay"  and  "prepare  variational  data  for  unit”.  In 
the  context  of  STAGSC-1  this  means  that  the  functions  performed  are 
preparatory  to  formulation  of  the  stiffness  matrix  during  execution 
with  STAGS2.  Figure  9.9  indicates  that  the  overlay  loads  the 
secondary  overlays  (7,1),  (7,2)  and  (7,3).  The  computational  flow  is 
controlled  mainly  by  subroutine  PREVU.  PREVU  processes  all  elements 
for  one  shell  unit  at  a  time.  It  calls  the  element  subroutines  GSBM, 
QUAF  and  TRINC  as  required.  These  subroutines  perform  all  necessary 
geometric  calculations,  strain-displacement  relationships, 
integration  point  coordinates,  weighting  factors  and  contributions  to 
the  mass  matrix. 

L0V7  -  loads  overlays  from  level  7 
0VE22  -  beam  element  overlay 

GSBM  -  master  routine  for  beam  element  generation 

MASSE  -  assembles  and  transforms  beam  element  mass  matrix 

MACUP  -  controls  formulation  and  update  of  element  constitutive  matrix 

0VE41  -  plate  element  overlay 

QUAF  -  master  routine  for  quadrilateral  plate  generation 
MAPXY  -  performs  bilinear  mapping 

FDF  -  finds  integration  weights  and  function  formulas  for  bilinear 
quadrilateral 

MSH  -  finds  coordinates  and  integration  points 
0VE31  -  triangular  element  overlay 

TRINC  -  master  routine  for  triangular  element  generation 

Figures  9.9  through  9.14  contain  full  details  of  the  subroutine 
linkages  in  this  overlay  and  the  secondary  overlays. 
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E.  Overlay  (10,0)  -  Program  OVSV  -  Data  Base  Preservation 

This  overlay  is  a  data  management  program  to  organize  the 
preprocessed  information  in  mass  storage  in  readiness  for  execution 
by  STAGS2.  The  subroutine  linkages  are  shown  in  Figure  9.15. 

F.  Overlay  (12,0)  -  Program  OVSI  -  Shell  Unit  Intersection 

Overlay  (12,0),  (Figure  9.16)  has  the  sole  purpose  of  checking  the 
interconnections  between  shell  units  for  consistency  and,  where 
inconsistencies  are  detected,  removing  them  if  possible. 

G.  Overlay  (13,0)  -  Program  OVSL  -  Shell  Unit  Loads 

Figure  9.17  shows  the  subroutine  calls  for  this  overlay.  Program 
OVSL  has  a  single  call  to  subroutine  LOADS  which  generates  a  loads 
file  for  a  given  shell  unit.  LOADS  generates  consistent  nodal  forces 
from  applied  loads  and  computes  additions  to  the  load  vector 
corresponding  to  imposed  displacements.  Output  of  the  load  file  and 
mass  file  is  controlled  by  subroutines  LOADOP  and  MASSOP. 

H.  Overlay  (25,0)  -  Program  OVEL  -  Element  Unit  Loads 

Program  OVEL  (Figure  9.18)  performs  functions  similar  to  OVSL  for  the 
element  unit.  The  controlling  routine  is  LOADE  with  the  force  vector 
being  determined  by  FORCEE.  FORCEE  does  not  currently  have  as  much 
capability  as  FORCES  in  OVSL  since  it  does  not  process  distributed 
loads. 

This  completes  the  description  of  the  major  functions  in  program  STAGS1 . 
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3.3  STAGS?  -  Execution  Phase 


The  structure  of  STAGS2  follows  the  same  standard  overlaying  concept  as 
STAGS1.  There  are  eight  primary  overlays  and  four  secondary.  Figure  3.3 
shows  the  overlay  links.  The  secondary  overlays  are  listed  with  their 
functions  in  Figure  3.4. 

STAGS2  functions  partly  through  subroutine  CONTRL  (Figure  9.19)  and  partly 
through  direct  calls  to  the  overlays.  The  direct  calls  are  for  the  purpose  of 
determining  storage  requirements  while  CONTRL  directs  all  analysis  functions. 
Also,  STAGS2  loads  user  subroutines  required  during  execution  (LUSER2). 

Oetails  of  the  subroutine  calls  for  STAGS2  are  contained  in  the  Appendix 
(Figures  9.19  through  9.39). 

A.  Overlay  (1,0)  -  Program  0V10  -  Data  Transfer  from  Pre-processor 

This  program  performs  the  many  complex  operations  required  to  start 
execution,  either  for  a  new  problem  or  a  restarted  problem.  Its  main 
functions  are  performed  by  three  subroutines;  ALL0C2,  DATAIN,  SETPAR, 
and  RSTRT.  ALL0C2  is  itself  a  complex  routine  which  determines  block 
sizes,  assigns  files,  sets  pointers  for  various  operations  (e.g., 
stiffness  matrix  computation),  determines  working  space  in  core  and 
initializes  file  manager  parameters.  Subroutine  DATAIN  initializes 
parameters.  SETPAR  also  performs  parameter  initialization,  e.g., 
initial  conditions  for  a  transient  analysis.  RSTRT  controls  a 
restart  analysis.  Subroutine  links  are  shown  in  Figures  9.22  and 
9.23. 


B.  Overlay  (2,0)  -  Program  0V20  -  Stiffness  Matrix  Decomposition 

Program  0V20  performs  one  of  the  crucial  stages  in  the  STAGSC-1 
analysis,  i.e.,  factorization  of  the  stiffness.  It  first  assembles 
the  total  stiffness  matrix  from  the  element  stiffness  file  in 
subroutine  ASEM.  ASEM  calls  ASEM2  which  adds  the  element 
contributions.  However,  the  actual  operations  are  carried  out  by 
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calling  a  COMPASS  coded  subroutine  ASEMA  (Figure  9.24)  which  renders 
the  details  somewhat  inaccessible  to  a  reviewer,  lhe  decomposition 
into  upper  and  lower  triangles  is  controlled  by  subroutine  FAG IUR 
with  the  actual  reduction  occurring  in  subroutine  FAGMI). 

C.  Overlay  (3,0)  -  Program  0V30  -  Eigenvalue  Solution 

The  eigenvalue  solver  is  probably  the  most  complex  program  in 
STAGSP.  Figure  9.26  provides  details  of  the  subroutine  links.  The 
major  functions  are  driven  by  subroutine  S1M1 1  which  controls  the 
computational  flow  for  simultaneous  iteration  for  a  cluster  of 
eigenvalues.  S1M1T  determines  the  number  of  eigenvectors  required  as 
a  subspace  for  the  simultaneous  inverse  iteration,  performs  the 
inverse  interation  and  solves  for  the  reduced  set  of  eigenvectors 
using  the  Householder  tridiagonal ization  and  QL  method.  The  major 
subroutines  in  which  these  computations  are  performed  are  EIGEN 
(Figure  9.26)  and  SOLVE.  EIGEN  calls  subroutines  TKE02  and  TQL2 
which  are  FORTRAN  versions  of  ALGOL  procedures  originally  developed 
by  Wilkinson,  Martin  and  Reinsch. 

0.  Overlay  (4,0)  Program  0V40  -  Formulation  of  Stiffness 

The  comments  provided  in  the  subroutines  called  by  this  overlay  refer 
to  "second  variation  of  strain  energy".  This  terminology  refers  to 
all  the  functions  normally  associated  with  the  generation  of 
stiffness  matrices  and  this  is  indeed  the  function  of  the  overlay. 
Figures  9.27,  9.28  and  9.29  show  the  subroutine  links  for  the 
overlay.  The  major  routines  called  are  VAR2,  CVR2,  VR2  and  VRbATA. 
VAR2  and  CVR2  are  the  controlling  routines  which  call  VRDATA  and 
VR2 .  VM DA T A  provides  all  the  necessary  information  specific  to  the 
element  type  while  VR2  performs  the  actual  calculations  of  the 
stiffness  contributions  at  a  given  integration  point.  Nonlinear 
terms  are  handled  by  VR2  and  also  effects  due  to  live  pressure 
loads.  Brief  functional  descriptions  of  significant  routines  are  as 
fol lows : 
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SDRM 
FDRV 

QUADF 
E32LL 

TRIDUV 

TRIDW2 

PENAL 

VR12D,  VR22D  -  computes  stiffness  terms  for  one  and  two-dimensional 
elements 

E.  Overlay  (5,0)  -  Program  0V50  -  First  Variation  of  Strain  Energy 

The  form  and  function  of  this  overlay  are  geared  directly  to  the 
method  of  solution  of  the  nonlinear  equations.  The  modified 
Newton-Raphson  method,  as  described  in  Section  4.5,  solves  the 
nonlinear  system  of  equations  iteratively  and  obtains  the  incremental 
displacement  vector  by  the  solution  of  the  equation 


-  Forms  strain-displacement  matrix  for  a  beam 

-  finds  interpolating  polynomials,  1st  and  2nd  partial 
derivatives  for  quadrilateral  elements 

-  performs  integrations  for  quadr i 1 ateral  s 

-  generates  shape  functions  for  live  pressure  loads  on 
triangular  elements 

-  finds  membrane  strain-displacement  matrix  for 
triangular  elements 

-  finds  curvature-displacement  matrix  for  triangular 
elements 

-  adds  penalty  terms  to  stiffness  matrix 


‘W  -  {V  =  -  [K(xm)rl  {K(xn)  <V  “  {R})  3J 

The  nonlinear  stiffness  K(x  )  is  determined  in  overlay  (4,0) 

m 

(second  variation  of  strain  energy).  This  overlay  forms  the  product 
K(xn)(xn)  directly  as  the  first  variation  of  strain  energy 
and  subtracts  the  contributions  from  the  force  vector  {R}.  The 
solution  algorithm  is  discussed  more  fully  in  Section  4.5. 

Figures  9.30,  9.31  and  9.32  show  the  subroutine  links  for  0V50.  The 
major  subroutines  called  are  ITER,  CVR1,  VR1  and  SOLVE.  ITER 
controls  the  interations  and  checks  convergence.  The  computation  of 
the  first  variation  is  controlled  by  CVR1.  This  pulls  in  all  the 
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necessary  element  information  (VkUATA)  required  for  computing  the 
first  variation  and  performs  ttie  computation  at  an  integration  point 
in  VK1.  Plasticity  calculations  are  performed  by  the  secondary 
overlay  (5,1)  in  program  CJV51.  CVkl  loads  program  UVbl  (figure 
9.31),  which  has  major  subroutines  as  follows: 

PLASTC  -  controls  plasticity  calculations  at  each  integration  point 
PLAST  -  performs  plastic  stress  calculations 

tquation  solving  is  performed  in  subroutines  SUL VI  and  SWEEP. 

F.  Overlay  (6,0)  -  Program  0V60  -  Dynamic  Response 

Program  0V60  performs  the  same  function  for  dynamic  response 
calculations  that  0V50  does  for  static  analysis.  The  procedures  are 
necessarily  more  complex  because  of  the  time  integration.  Figures 
9.33  and  9.34  show  the  subroutine  links.  The  controlling  subroutine 
for  0V60  is  UYNR  which  calls  CVR1,  ODES  and  EXPLC.  CVR1  is  the  same 
as  0V50  and  performs  the  same  function,  i.e.,  sets  up  the  right,  hand 
side  vector  for  the  nonlinear  solution  as  the  first  variation  of 
strain  energy.  The  dynamic  response  obtained  using  the  explicit 
integration  operator  does  not  require  the  solution  of  equations  anu 
is  therefore  called  directly  by  DYNR.  Solutions  using  implicit 
integration  operators  are  performed  by  subroutine  DUES.  DDLS  is  a 
general,  multi-step  ordinary  differential  equation  solver  and  is 
called  once  per  time  step.  It  controls  the  time  step  either  in  the 
automatic  mode  or  in  the  "fixed"  mode  in  which  it  only  intervenes  it 
the  time  step  needs  to  be  decreased.  It  also  provides  starting 
procedures  and  handles  damping.  ODES  calls  MSTEP,  SOLVE,  NITER  and 
NEXT.  MSTEP  computes  the  predictor-corrector  formulae  for  the  four 
implicit  schemes.  There  is  also  a  link  to  a  user  defined  operator  by 
a  call  to  a  user  written  subroutine  USTEP.  This  capability  is  not 
currently  documented  and  is  not  included  in  LUSER3.  NITER  controls 
the  iterations  for  nonlinear  equation  solving  with  SOLVE  providing 
the  solution  procedure.  Subroutine  NEXT  provides  the  coding  for 
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time  step  control.  0V60  calls  a  separate  overlay  for  plasticity, 
OVlbl  in  dynamic  calculations,  but  this  appears  to  be  identical  with 
the  plasticity  routine  used  in  static  analysis  0VL51  (Figure  9.3 2). 

G.  Overlay  (7,0)  -  Program  0V70  -  Solution  Strategy. 

Program  0V70  is  called  by  subroutine  CONTRL  to  control  the  static 
nonlinear  solution  strategy.  The  controlling  subroutine  within  0V70 
is  DATA1  which  controls  the  output  for  each  load  step  and  calls  the 
major  subroutines  STRAT,  EQCHK,  SOATA  and  OUTSLLl.  STRAT  controls 
load  step  size  and  adjusts  it  if  necessary  depending  on  the  rate  of 
convergence.  It  also  controls  refactoring  of  the  stiffness  and 
extrapolation  of  displacements  for  the  next  load  step.  Subroutine 
EQCHK  performs  an  overall  equilibrium  check  but  is  disabled  in  the 
version  evaluated  since  the  programming  is  not  yet  complete.  SDATA 
maintains  the  solution  data  file  (SOD),  plastic  stress  history  and 
also  writes  TAPE8  for  nonlinear  analysis  in  which  periodic 
eigen-solutions  may  be  derived.  Thus,  estimates  of  bifurcation 
buckling  loads  may  be  obtained  at  various  points  in  a  nonlinear 
loading  history.  Finally,  subroutine  OUTSLD  controls  output  of 
selected  displacements. 

H.  Overlay  (10,0)  -  Program  0V80  -  Stress  Computation  and  Output 

Overall  control  of  stress  and  strain  computation  and  output  is 
exercised  through  this  overlay.  A  master  routine  SIGMA  calls 
secondary  overlays  0V81  and  0V82  and  major  subroutines  OUTSLS,  SRES, 
VRDATA  and  PREFAB.  SIGMA  controls  stress  calculations  element  by 
element.  Element  data  is  brought  in  by  means  of  VRDATA.  SRES 
controls  the  actual  calculation  of  strains  and  stress  resultants. 

The  secondary  overlays  0V81  and  0V82  do  the  element-specific 
computations;  0V81  handles  1-D  (beam)  elements  while  0V82  deals  with 
the  2-D  elements.  Figures  9.37,  9.38  and  9.39  provide  the  details  of 
the  subroutine  1  inks. 
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3.4  DATA  MANAGEMENT 


The  storaqp  and  retrieval  of  data  in  STAGSC-1  is  accomplished  by  two  separate 
schemes,  t  irst  ,  there  is  tin*  need  to  transfer  large  files  for  vector  and 
matrix  operations  between  core  and  mass  storage.  The  size  of  the  files  may  be 
particularly  large  when  time  integration  is  being  performed  and  vectors  must 
be  made  available  for  several  time  steps.  Also,  during  an  eigensolution, 
twenty  or  thirty  subspace  vectors  are  being  manipulated  together  with  mass  and 
stiffness  matrices.  Second,  there  is  the  need  to  obtain  relatively  small 
amounts  of  data  from  tables  in  order  to  generate  element  stiffness  matrices  or 
calculate  element  stresses.  These  tables  are  themselves  lengthy  and  the  mode 
of  retrieval  may  be  described  as  quasi-random  access;  this  is  because  transfer 
of  successive  sections  of  data  may  be  from  regions  of  the  file  which  are 
adjacent . 

3.4.1  FILE  MANAGER  -  FMM 

FMM  is  a  routine  designed  to  manage  working  space  in  blank  common  for  vector 
and  matrix  manipulations.  When  FMM  is  called,  the  argument  first  identifies 
the  number  of  files  to  be  located  simultaneously  in  core,  a  list  of  the  file 
numbers  and  their  lengths  and  also  a  priority  indication  which  says  whether 
the  file  is  to  be  saved  or  not  (in  mass  storage)  when  removed  from  core.  FMM 
provides  as  output  the  address  in  blank  common  of  each  file.  During 
execution,  FMM  checks  if  a  given  file  is  already  in  core;  determines  whether 
the  space  required  by  the  file  is  available;  adds  the  file  to  core  and 
performs  other  housekeeping  functions.  The  present  version  does  not  utilize 
the  LCM  (large  core  memory)  feature  available  on  the  CDC  7600. 

FMM  also  makes  use  of  a  number  of  other  utility  routines  for  performing 
specific  operations.  The  subroutine  calls  may  be  found  in  the  Appendix  in 
F igure  9.3. 
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3.4.^  VIRTUAL  MEMORY  OPERATIONS  -  10GEIR. 


The  transfer  of  tabular  data  from  mass  storage  into  core  is  accomplished  using 
a  buffering  technique.  Such  data  are  used  in  element  matrix  generation, 
stress  calculations,  etc.,  and  are  stored  as  tables  in  one  lengthy  file.  On 
the  other  hand,  the  data  are  needed  only  in  relatively  small  blocks  at  any 
time.  The  technique  used  in  writing  the  file  is  to  divide  it  into  a  number  of 
records  of  convenient  length.  The  record  length  is  chosen  so  that  typically  6 
to  10  records  can  be  accomodated  in  core  at  a  given  time.  This  process  is 
reasonably  efficient  since  the  required  data  are  often  in  adjacent  blocks  if 
not  all  in  one  block.  The  controlling  subroutine  is  I06ETR,  which  searches 
the  buffer  for  the  required  record  and  reads  it  in  if  it  is  not  already  there 
(having  checked  for  space  availability).  If  the  buffer  is  full,  the  last 
record  is  evicted  and  written  to  mass  storage.  IOGETR  calls  a  few  utility 
subroutines  and  the  links  are  shown  in  Figure  9.3. 
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Figure  3.2  Secondary  Overlay  Structure  of  STAGS  1 
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4.0  FUNCTIONAL  DESCRIPTION 


The  purpose  of  this  section  is  to  provide  a  detailed  overview  of  the  major 
capabilities  and  analytical  methods  used  in  STAGSC-1. 

4.1  ANALYSIS  OPTIONS 

As  outlined  in  the  Introduction  (Section  1),  STAGSC-1  is  a  general  purpose, 
thin-shell,  structural  analysis  program,  designed  principally  for  the 
nonlinear  static  and  dynamic  analysis  of  thin  shells.  There  are  seven 
different  analysis  options  available  to  the  user. 

o  Linear  static  analysis 

o  Bifurcation  buckling  analysis  from  a  linear  stress  state 
o  Small  vibrations  (stress  free  state) 

o  Nonlinear  static  analysis 

o  Bifurcation  buckling  analysis  (nonlinear  stress  state) 
o  Small  vibrations  (linear  or  nonlinear  stress  state) 
o  Transient  response  analysis  (linear  or  nonlinear) 

4.1.1  LINEAR  STATIC  ANALYSIS 

Although  there  are  numerous  general  purpose  finite  element  programs  which 
provide  thin  shell  elements  for  linear  static  analysis,  STAGSC-1  has  a  number 
of  features  which  make  it  an  attractive  choice  for  this  application. 
Specifically,  these  are:  (i)  built-in  geometries  for  regular  shell  surfaces 
such  as  cylinder,  cone,  flat  plate,  torus,  sphere,  etc.;  (ii)  user  subroutine 
capability  for  defining  surface  deviation  with  respect  to  some  reference 
surface;  (iii)  discrete  surface  stiffeners  (both  orthogonal  and  skewed);  and 
(iv)  multilayer  shell  wall  construction.  Thus,  complex  shell  geometries  can 
be  modeled  with  relative  ease. 
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Loads  may  be  applied  directly  at  mesh  points  (defined  by  “row"  and  " column " 
»umbers--Sec t ion  4.2)  as  forces  and/or  moments,  as  line  loads  or  moments  or  as 
surface  tractions.  Their  directions  may  be  specified  in  global  or  local 
surface  directions.  Thermal  loads  can  also  be  generated  by  specifying 
reference  surface  temperatures  in  a  user  coded  subroutine.  Temperature 
variation  in  a  stiffener  cross-section  may  be  prescribed  (but  no  variation  is 
permitted  through  the  shell  wall  thickness).  Body  forces  can  be  specified  by 
means  of  acceleration  vectors  in  both  translation  and  rotation.  Displacement 
boundary  conditions  may  be  applied  as  discrete  constraints  at  interior  or 
boundary  mesh  points  or  by  specialized  conditions  (e.g.,  simple  supports, 
clamped,  etc.)  along  shell  boundary  edges. 

The  range  of  basic  capabilities  for  static  analysis  is  therefore  quite 
adequate  for  simple  stress/displacement  analysis  of  thin  shells  and  includes 
unique  features  (such  as  the  geometric  deviations  from  a  reference  surface) 
which  are  a  definite  incentive  for  its  practical  use. 

4.1.2  BIFURCATION  BUCKLING  ANALYSIS 

The  analysis  of  bifurcation  buckling  of  shells  has  been  an  important 
capability  in  STAGS  since  it  was  introduced  in  an  early  version  in  about  1970 
(see  Figure  1.1).  This  part  of  the  program  is  therefore  one  of  the  most 
hiqhly  developed  and  would  probably  be  a  common  reason  for  selection  by  a 
potential  user. 

Bifurcation  buckling  can  be  investigated  for  structures  which  have  linear  or 
nonlinear  prebuckling  stress  states.  An  example  of  linear  prestress  is  given 
by  a  flat  plate  with  in-plane  loading  only.  Bifurcation  buckling  is  then 
defined  as  the  value  of  the  load  at  which  a  laterally  displaced  configuration 
can  also  be  in  equilibrium  (secondary  loading  path).  A  shell  of  revolution, 
such  as  a  shallow  spherical  cap,  will  exhibit  nonlinear  prebuckling  behavior. 
Under  antisymmetric  loading,  the  shell  softens  or  stiffens  depending  on  the 
loadinq  direction  (e.g.,  external  or  internal  pressure).  Bifurcation  buckling 
may  then  occur  as  a  nonsymmetric  deformation  mode.  Figure  4.1  shows  the 
characteristic  behavior  of  systems  which  exhibit  linear  and  nonlinear 
prebuckling  primary  load  paths. 
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For  a  simple  linear  prebuckling  problem,  the  eigenvalue  analysis  provides  the 
multiplier  by  which  the  basic  applied  load  system  must  be  factored  to  obtain 
the  critical  load  level.  STAGS  allows  for  the  specification  of  two 
independent  load  systems  {F^l  and  (Fg)  which  are  characterized  by 
multipliers  and  Pg.  Thus,  for  example,  a  cylinder  with  internal 
pressure  and  axial  compression  will  require  the  axial  load  to  be  designated  as 
System  A  and  the  pressure  as  System  B.  For  a  given  pressure  loading,  the 
total  magnitude  of  the  two  load  systems  will  be 

{F)T0TAL  =  *  PA  lFA}  +  PB  {FB}  4,1 

where  a  is  the  eigenvalue  at  bifurcation. 

The  case  of  bifurcation  buckling  analysis  with  a  nonlinear  prebuckling  stress 
state  is  an  option  with  far  less  general  applicability.  However,  the  more 
sophisticated  user  who  needs  to  perform  a  nonlinear  collapse  analysis  for  a 
general  shell  can  take  advantage  of  a  number  of  subtleties  which  the  nonlinear 
bifurcation  capability  provides.  The  STAGS  theoretical  manual  [3]  offers  a 
very  detailed  and  thorough  discussion  of  bifurcation  buckling  and  collapse 
analysis.  Almroth  and  Brogan  [5]  give  a  number  of  examples  in  which  the 
nonlinear  collapse  loads  are  calculated  and  compared  with  bifurcation  buckling 
loads  obtained  using  linear  prebuckling  analysis.  It  is  shown  the  linear 
bifurcation  loads  may  be  greater  or  less  than  actual  nonlinear  collapse 
loads.  For  example,  an  elliptic  cone  undergoing  uniform  end  shortening  will 
collapse  at  a  load  over  twice  that  predicted  by  linear  bifurcation,  while  a 
cylindrical  panel  with  its  ends  simply  supported  ("Venetian  blind"  model)  will 
collapse  at  a  load  five  times  smaller  than  the  bifurcation  load.  The 
usefulness  of  the  nonlinear  bifurcation  analysis  appears  to  be  in  its 
application  as  an  adjunct  to  a  full  nonlinear  collapse  analysis.  If  nonlinear 
collapse  analysis  is  performed  on  an  imperfection  sensitive  structure,  the 
analysis  may  fail  due  to  ill-conditioning  of  the  equations  at  some  load  step. 

A  bifurcation  buckling  analysis  carried  out  at  this  point  will  yield  a 
buckling  mode  which  will  indicate  the  type  of  imperfection  which  will  direct 
the  solution  into  the  secondary  path  (see  Figure  4.1(b)).  Another  practical 


38 

09I8B-84B:? 

(S3034)  3 


situation  is  where  a  nonlinear  analysis  is  carried  out  only  up  to  a  design 
load.  A  nonlinear  bifurcation  analysis  at  this  point  will  provide  an  estimate 
of  the  remaining  margin  to  collapse. 


4.1.3  SMALL  VIBRATION  ANALYSIS 

The  capability  of  STA6SC-1  with  respect  to  vibration  mode  analysis  is  very 
similar  to  the  bifurcation  buckling  capability.  An  eigenvalue  solution  for 
vibration  modes  and  frequencies  may  be  obtained  for  a  stress-free  structure  or 
for  a  linear  or  nonlinear  stress  state.  This  is  particularly  relevant  for  the 
analysis  of  shell  structures  where  the  presence  of  pressure  loading  (internal 
or  external  to  the  shell)  is  common  prior  to  dynamic  loading. 

In  the  case  of  vibration  all  eigenvalues  will  be  positive,  whereas  in  the  case 
of  bifurcation  buckling  negative  eigenvalues  may  be  obtained  (e.g.,  in  a  shear 
loaded  plate) . 

In  both  the  bifurcation  and  vibration  analysis  options,  it  is  possible  to 
define  by  input  a  uniform  stress  state  directly  in  terms  of  direct  and  shear 
stress  resultants  as  an  alternative  to  the  generation  of  such  stress  states  by 
means  of  applied  loads,  displacements  or  temperatures. 

STAGSC-1  also  permits  the  user  to  specify  concentrated  (lumped)  masses 
directly  at  node  points.  It  should  be  noted  that  lumped  rotational  inertias 
cannot  be  specified. 

A  further  limitation  appears  to  be  that  vibration  modes  for  unsupported 
structures  ("free-free")  cannot  be  determined  because  the  eigensolver  needs  to 
solve  the  system  equations.  Rigid  body  constraints  will  permit  a  solution  to 
be  obtained. 

4.1.4  NONLINEAR  STATIC  ANALYSIS 

Both  geometric  (large  displacements)  and  material  nonlinearities  can  be 
included  in  a  STAGSC-1  analysis.  The  geometrically  nonlinear  analysis  is 
based  on  the  following; 
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A.  Small  strain  based  on  engineering  stress  and  strain  relationships 

(c  <0.1) 

B.  Strain-displacement  equations  retaining  nonlinear  rotation  terms 
(moderate  rotations  <0.3  radians) 

C.  Incremental  solution  of  equations  with  iterations  within  each 
increment  using  the  modified  Newton-Raphson  method. 

Material  nonl inearities  considered  are  due  to  plasticity  only.  No  creep  or 
viscoelastic  behavior  is  incorporated  in  STAGSC-1.  Plasticity  is  handled 
according  to  the  White-Bessel ing  theory  [7].  This  is  equivalent  to 
elastic-perfectly  plastic  behavior,  bilinear  kinematic  hardening  or 
multilinear  hardening  depending  on  the  number  of  plastic  parameters 
specified.  The  theory  is  outlined  in  greater  detail  in  Section  4.4. 

The  plastic  strains  are  computed  and  used  to  generate  pseudo-force  vectors, 
i.e.,  an  initial  strain  method  is  implemented. 

4.1.5  TRANSIENT  RESPONSE  ANALYSIS 

This  capability  represents  one  of  the  major  strengths  of  the  STAGSC-1  program 
since,  like  the  bifurcation  analysis,  it  has  been  under  development  for  a 
number  of  years.  The  program  can  solve  transient  problems  with  a  wide  range 
of  excitation  using  one  of  five  different  transient  integration  operators. 
System  damping  may  be  introduced  as  Rayleigh  viscous  damping,  with  constant 
stiffness  and  mass  matrix  multipliers,  plus  an  additional  contribution  from 
velocity  dependent  forces.  In  addition,  the  full  range  of  nonlinearities 
available  for  static  analysis  can  be  utilized  in  transient  response. 

Forcing  functions  may  be  specified  in  terms  of  nodal  loading  or  displacement 
patterns,  with  time  dependencies  either  (a)  according  to  certain  specified 
formats  or  (b)  input  through  a  user  coded  subroutine  (FORCET).  The  specified 
formats  are: 
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o  piecewise  linear  function 
o  tr iqonometric  function 
o  exponential  decay 

In  addition,  initial  velocities  and  displacements  can  be  specified  for  any 
degrees  of  freedom. 

The  integration  operators  implemented  in  STAGSC-1  are  as  follows: 
o  explicit  (central  difference) 

o  implicit  (trapezoidal.  Gear  2nd  and  3rd  order,  and  Park's  method) 

The  time  step  for  the  central  difference  method  is,  of  course,  fixed  and  must 
he  selected  by  the  user.  For  the  implicit  methods,  either  a  fixed  time  step 
may  be  used  or  there  is  an  internal  algorithm  for  automatic  time  step 
control.  The  automatic  feature  is,  however,  presently  regarded  as 
experimental  by  the  developers. 

4.2  SURFACE  AND  MESH  GEOMETRY 

The  STAGSC-1  philosophy  for  developing  the  shell  geometry  and  finite  element 
discretization  contains  some  rather  unfamiliar  concepts  and  terminology.  A 
central  notion  is  that  of  the  so-called  "shell  unit".  This  can  refer  to  the 
description  of  a  specific  portion  of  the  shell  surface  or  to  the  whole 
surface.  The  complete  shell  may  be  defined  using  up  to  thirty  shell  units. 

In  addition,  or  as  an  alternative,  the  structure  can  be  defined  in  terms  of  an 
element  unit.  The  basic  distinction  between  these  two  concepts  is  that  the 
shell  unit  defines  a  geometric  surface  with  a  rectangular  grid  work  mapped 
onto  the  surface,  while  the  element  unit  is  actually  a  direct  assemblage  of 
elements  which  may  or  may  not  define  a  shell. 

The  shell  unit  grid  is  then  overlaid  with  a  mesh  of  finite  elements  which  may 
use  all  qridpoints  of  the  shell  unit  or  only  a  subset.  This  concept  allows 
the  use  of  a  library  of  standard  shell  geometries  (e.g.,  cylinder,  sphere. 
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torus,  etc.)  which  may  then  be  adapted  to  model  cutouts  or  stiffened  areas  by 
the  appropriate  omission  or  addition  of  elements.  Discrete  stiffeners  may 
also  be  attached  along  arbitrary  paths  on  the  shell  surface. 

4.2.1  SURFACE  GEOMETRY  LIBRARY 

There  are  eleven  standard  geometries  in  STAGSC-1  which  are  selected  using  the 
input  key  ISHELL.  The  geometries  are  listed  in  Table  4.1 

Each  geometry  has  four  edges,  two  of  which  may  be  subsequently  joined  to  each 
other  to  form  a  closed  surface  (except  for  the  rectangle  and  quadrilateral). 
Alternatively,  edges  may  be  joined  to  other  shell  units  or  may  have  boundary 
conditions  applied. 

4.2.2  USER  DEFINED  SURFACE  GEOMETRY 

This  may  be  accomplished  using  the  user-coded  subroutine  LAME.  This 
subroutine  may  define  global  coordinates  for  the  surface  and  their  first  order 
derivatives  if  flat  elements  are  being  used.  This  is  currently  the  only 
usable  option  but,  anticipating  the  introduction  of  curved  elements,  the  user 
can  define  directly  the  coefficients  of  the  first  and  second  fundamental  forms 
or,  alternatively,  all  the  derivatives  necessary  for  internal  computation  of 
the  coefficients. 

4.2.3  SURFACE  GRID  AND  ELEMENT  MESH 

Given  that  a  reference  surface  geometry  has  been  defined,  a  gridwork  must  be 
mapped  onto  the  surface.  In  its  simplest  form,  this  is  accomplished  by 
specifying  numbers  of  rows  and  columns  which  generates  a  regular  gridwork  in 
terms  of  the  surface  coordinates.  Additional  options  are: 

o  irregular  grid  spacing  by  means  of  definition  of  different  segments 
o  grid  definition  by  user  subroutine  IUGRID 

The  element  mesh  is  conceptually  distinct  from  the  surface  grid.  This  is 
specified  separately  and  may  be  defined  in  a  number  of  ways.  In  its  simplest 
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form,  the  element  mesh  is  identical  with  the  surface  fjrid  in  that  all  grid 
points  define  element  nodes.  Alternatives  are: 

A.  Irregular  mesh  in  which  cutouts  can  be  defined. 

B.  Specialized  mesh  in  which  grid  points  are  used  selectively  and  a  mesh 
of  varying  refinement  can  be  obtained. 

C.  Subregion  or  "patch"  concept  in  which  groups  of  elements  are  defined 
within  certain  row  and  column  boundaries.  This  allows  use  of 
different  element  types  in  different  regions  of  the  shell  surface. 

4.2.4  ELEMENT  UNITS 

An  element  unit  can  be  defined  as  a  “stand-alone"  unit  which  describes  the 
total  structure  or  it  may  be  used  in  conjunction  with  a  shell  unit.  In  the 
latter  case,  element  unit  nodes  can  be  nodes  on  the  shell  unit  or  separately 
defined  (auxiliary)  nodes  or  a  combination.  Thus,  parts  of  the  structure 
which  are  not  shell-like  can  be  connected  to  the  shell.  This  is  in  addition 
to  the  capability  which  exists  in  shell  units  for  specifying  stiffeners  on  the 
shell  surface. 

Node  geometry  for  the  element  unit  must  be  defined  individually  for  each  node 
or  by  means  of  the  user  subroutine  USRPT  or  by  a  combination  of  both.  No 
other  options  are  available. 

The  directions  of  the  degrees  of  freedom  can  be  separately  defined  for  the 
auxiliary  nodes  in  the  element  unit. 

4.3  FINITE  ELEMENT  LIBRARY 

The  core  of  the  finite  element  library  in  STAGSC-1  is  the  series  of  triangular 
and  quadrilateral  shell  elements  based  on  the  Clough-Fel ippa  quadrilateral 
bending  element  [8].  Since  these  are  all  flat  elements,  the  actual  curved 
shell  geometry  is  always  approximated  by  a  faceted  surface.  This  has 
implications  for  interelement  compatibility  which  will  be  discussed  later. 
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In  addition  to  the  shell  elements  there  are  membrane  versions  of  the  elements 
and  also  a  series  of  beam  elements  designed  to  be  compatible  with  the  shells. 
STAGSC-1  also  has  provision  for  a  series  of  general  linear  and  nonlinear 
springs  and  also  transition  membrane  and  plate  elements.  The  latter  are 
designed  for  regions  of  changing  mesh  refinement.  Neither  the  springs  nor  the 
transition  elements  were  implemented  in  the  version  of  the  program  under 
evaluation.  In  total,  there  are  twenty  elements  currently  available  in 
STAGSC-1,  of  which  six  are  beam  elements  and  fourteen  are  membrane  and  shell 
elements.  Table  4.?  provides  a  summary  description  of  the  elements 
implemented. 

The  bar  and  beam  elements  are  fairly  standard  except  for  220  and  221  which 
have  quadratic  shape  functions  for  twist.  The  inclusion  of  the  center  node 
makes  them  compatible  with  the  majority  of  the  shell  elements  in  STAGSC-1. 

Also,  the  center  node  provides  better  results  when  displacements  (rotations) 
are  relatively  large.  For  a  general  thin  shell  analysis  program  such  as 
STAGSC-1,  the  single  most  important  aspect  of  the  program  has  to  be  the 
properties  and  performance  of  the  shell  elements  themselves.  As  has  already 
been  mentioned,  STAGSC-1  does  not,  at  present,  have  available  a  curved  shell 
element  and  this  introduces  inevitable  incompatibilities.  However, 
considerable  effort  has  been  spent  by  the  developers  on  minimizing  these 
shortcomings . 

The  basis  for  both  the  triangular  and  quadrilateral  elements  is  the 
Clough-Fel ippa  triangle  (LCCT-12  in  Reference  8).  This  is  a  triangular 
bending  element  consisting  of  3  sub-elements  with  interior  nodes  condensed 
out.  The  lateral  displacements  have  therefore  a  piecewise  cubic 
distribution.  The  addition  of  in-plane  degrees  of  freedom  and  membrane  shape 
functions  gives  rise  to  the  quadrilateral  element.  Reference  8  describes  a 
quadrilateral  bending  element  (0-19)  which  is  derived  from  four  LCCT-12 
elements  with  the  internal  freedoms  condensed  out  and  the  mid-side  rotations 
on  the  four  other  edges  constrained  to  be  the  average  of  the  adjacent  nodal 
components.  The  420  (QUARC)  elements  are  variants  of  the  Q-19  element  with 
the  addition  of  translat ional  freedoms  at  the  mid-side  nodes  parallel  and 
normal  to  the  edges.  Table  4.2  gives  details  of  the  resulting  shape  functions. 
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The  410  (QUAF)  elements  were  specially  developed  to  remove  the  displacement 
incompatibility  which  exists  when  flat  shell  elements  do  not  lie  in  one 
plane.  If  edge  displacements  are  to  be  compatible,  then  the  transverse 
displacement  shape  functions  must  be  of  the  same  order  as  the  in-plane 
functions.  The  in-plane  functions  must  therefore  be  cubic.  This  was 
accomplished  in  element  410  by  introducing  a  normal  rotation  at  each  of  four 
corners.  Element  411  carried  this  one  step  further  by  introducing  2  rotations 
at  each  corner,  thus  permitting  individual  rotation  of  each  adjacent  side  and 
thereby  permitting  shear  strain  at  the  corner.  In  addition,  tangential 
displacements  at  mid-side  nodes  are  also  included  in  this  element. 

4.4  CONSTITUTIVE  RELATIONSHIPS 

The  number  of  constitutive  behavior  models  available  in  STAGSC-1  is  somewhat 
limited.  Orthotropic  elastic  behavior  and  plasticity  are  the  two  major 
options.  Creep  or  viscoelasticity  are  not  available  in  the  present  version  of 
STAGSC-1  .* 

Elastic  properties  are  specified  with  respect  to  principal  material  directions 
for  an  orthotropic  material.  Specialization  to  the  isotropic  case  is 
trivial.  Plasticity  is  based  on  the  White-Bessel ing  (mechanical  sublayer) 
model.  The  theoretical  manual  [3]  discusses  various  types  of  plastic 
constitutive  behavior  and  describes  the  White-Bessel ing  model  in  some  detail. 
An  advantage  of  the  W-B  model  is  that  the  uniaxial  stress-strain  curve  can  be 
represented  with  fair  accuracy  by  choosing  a  sufficient  number  of  components 
(or  "sublayers") .  The  minimum  number  of  components  (2)  automatical ly  yields 
the  bilinear  kinematic  hardening  theory. 

The  input  of  material  properties  is  based  on  the  specification  of  a  material 
type  number.  Up  to  30  different  material  types  may  be  specified.  In  this  way 
different  material  properties  can  be  assigned  to  different  regions  of  the 
structure  (as  defined  by  different  shell  units).  Also,  mutilayered. 


♦Development  of  a  creep  version  is  being  sponsored  by  NSRDC.  The  creep 
implementation  is  to  be  the  same  as  in  the  BOSOR  4  program  [4]. 
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composite-shel 1-wal 1  construction  can  be  simulated  by  defining  different 
orthotropic  properties  for  each  layer.  There  is  no  provision  for 
incorporating  temperature  dependence  of  material  properties.  Although  the 
users  manual  [2]  states  otherwise,  material  properties  cannot  vary 
continuously  through  the  wall  thickness. 

The  elastic  properties  may  also  be  input  through  a  user  written  subroutine 
WALL.  This  gives  the  user  the  ability  to  vary  the  elastic  coefficients 
continuously  over  the  shell  surface.  The  elastic-plastic  stress-strain  data, 
however,  can  only  be  input  via  cards.  Any  variation  throughout  the  structure 
must  be  defined  by  varying  the  material  type. 

In  summary,  the  constitutive  capability  in  STA6SC-1  is  adequate  for  shell 
structures  operating  in  an  environment  where  temperatures  are  below  the  creep 
range  for  the  material.  A  simplified  approximation  to  elastic  temperature 
dependence  could  conceivably  be  achieved  through  the  user  subroutine  WALL  by 
correlating  temperature,  spatial  position  and  elastic  properties. 

4.5  LINEAR  AND  NONLINEAR  ANALYSIS 

STAGSC-1  is  primarily  a  tool  for  nonlinear  analysis  although  a  purely  linear 
analysis  option  is  available  and  is  very  economical.  Therefore,  linear  static 
and  dynamic  analyses  may  be  appropriately  performed  for  shell  problems  using 
STAGSC-1  because  of  its  many  features  which  facilitate  the  analysis  of 
stiffened  or  unstiffened  shells  (see  Sections  4.1,  4.2,  4.9  and  4.13). 

However,  the  bulk  of  the  developmental  effort  behind  STAGSC-1  has  been  devoted 
to  its  nonlinear  solution  algorithms  and  to  the  implementation  of  a  finite 
element  library.  It  is  in  this  context  that  the  program  must  mainly  be 
discussed. 

The  basic  method  of  equation  solving  is  the  modified  Newton-Raphson  method 
(MNR),  with  periodic  updating  of  the  stiffness  ("refactoring"  in  STAGSC-1 
terminology).  The  program  is  self-adaptive  to  a  certain  extent  in  that  it  can 
switch  to  a  full  Newton-Raphson  method  (FNR)  if  indicated  by  convergence 
behavior.  It  must  be  pointed  out,  however,  that  the  handlinq  of 
nonlinearities  due  to  plasticity  is  by  the  initial  strain  method.  Thus,  the 
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stiffness  matrix  is  never  modified  to  reflect  changes  due  to  the  current 
material  state.  In  these  circumstances,  the  method  must  be  viewed  as  a  hybrid 
technique. 


In  order  to  discuss  the  solution  procedure  more  specifically,  some  basic 
equations  will  be  developed.  For  a  static,  nonlinear  problem  the  equations  of 
equilibrium  may  be  written  as  a  Taylor  expansion  of  the  total  force  vector 
[F ( x )  ]  (sum  of  applied,  restoring  and  residual  force  vectors)  about  the 
currently  deformed  state; 


{F(x  , . ) }  —  {F (x  )}  +  -~T  I  ,  ,  ,  ,  (  {x  ,,  }  -  (x  })  +  terms  of  higher 

n+1  n"  3{xl  (x>  =  {xni  n+1'  n"  or(jer 

=  lF^xn^}  +  HR  1  (x )=  {x  }  {AX)  +  terms  of  hi9her  order 
=  0  4.2 

3(F) 

In  this  notation,  the  derivative  ■■■  '  |  is  the  negative  of  the 

B  )  {X  )-  {X^  } 

nonlinear  stiffness  matrix  [K(xn)]. 

A  fundamental  concept  of  the  STAGSC-1  program  is  to  treat  the  product  vector 
[K(x)]ux)  as  a  nonlinear  operator  L  acting  on  the  incremental 
displacements.  Thus,  equation  4.2  becomes: 


tF ( xn  +  i )  >  =  {L ( x ^ ) >  -  {R }  +  (hiqher  order  terms)  =  0  4.3 

where  the  operator  L  is  the  first  derivative  of  the  strain  energy  functional. 
The  nonlinear  stiffness  matrix  [K(x)]  is  then  the  first  derivative  of  L, 
[L'(x)].  The  significance  of  this  goes  far  beyond  the  formal  statement  of  the 
equilibrium  conditions  because  the  nonlinear  solution  algorithm  is  based  on 
the  direct  formulation  of  { L * ( x )  > .  This  has  important  advantages  in  the 
implementation  of  the  Newton-Raphson  method.  This  may  be  stated  in  the 
following  terms 


txn+l‘  =  lV  ‘  tF(*n5)  or 

ixn+i>  -  (xn>  =  [L'(xn)]_1  ( {R }  -  lL ( x n )  })  4.4 
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Fquation  4.4  implies  a  full  update  of  the  stiffness  matrix  and  re-solution 
(refactoring  in  STAGS  terminology)  at  each  iteration.  Basically,  a  modified 

I 

Newton-Raphson  procedure  is  implemented  which  performs  a  refactoring  of  the 
stiffness  only  when  indicated  by  convergence  criteria.  The  MNR  algorithm 
implemented  in  STAGSC-1  can  be  written  then  as 

Un+1}  '  {V  =  [L,(xm^_1  ( {R }  "  <L(xn) >)  4-5 

where  [L(x  1  is  represented  by  the  factored  stiffness  matrix  obtained  at 
some  previous  iteration  or  step.  The  operator  L(xn)  is,  however,  defined 
for  the  current  solution  vector  t xn i  since  it  is  formed  directly  as  a 
vector. 

The  nonlinearities  included  in  the  computation  of  L(xn)  are  purely  geometric 
and  plasticity  effects  are  handled  separately  as  pseudo-force  contributions  to 
the  loading  vector  tR > .  Thus,  plasticity  corrections  are  computed  after 
each  MNR  iteration  and  the  modified  vector  {R}  is  then  used  in  the  next 
iteration. 

Figure  4.2  illustrates  the  MNR  algorithm  and  plasticity  solution  in  an  overall 
sense  as  implemented  in  STAGSC-1.  It  should  be  pointed  out  that  this  is  a 
conceptual  flow  chart  and  does  not  represent  the  actual  program  flow. 

The  parameters  over  which  the  user  has  control  at  the  time  of  input  are 

o  total  number  of  times  the  step  size  may  be  cut  (NCUT) 
o  total  number  of  refactorings  allowed  (NEWT) 
o  initial  solution  estimate  (NSTRAT) 
o  convergence  tolerance  (DELEX) 
o  relaxation  factor  (WUND) 

By  the  use  of  a  negative  value  of  NEWT,  refactoring  can  be  enforced  at  desired 
load  step  intervals  (including  every  iteration).  Note:  the  program 
automatically  doubles  the  load  step  after  two  successive  steps  with  single 
iteration  convergence. 
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Tillerson,  Stricklin  and  Haisler  [9]  in  their  excellent  survey  of  numerical 
methods  for  solving  nonlinear  problems,  conclude  that  the  MNR  method  is 
probably  one  of  the  most  widely  used  for  general  nonlinear  problems  and  that 
for  structural  problems  it  is  best  suited  for  those  in  which  geometric 
nonlinearities  predominate.  Some  doubts  are  expressed  about  including 
material  nonlinearities  because  of  the  problem  of  unloading.  The  situation 
can  arise  using  MNR  where  elastic  unloading  is  not  correctly  handled  because 
the  factored  stiffness  corresponds  to  a  "tangent"  stiffness  based  on  a  prior 
loading  step.  On  the  other  hand,  Bushnell  [10]  describes  a  subincremental 
plasticity  fomulation  in  which  plasticity  calculations  are  performed  outside 
an  inner  Newton-Raphson  loop.  The  question  of  plastic  unloading  was  not 
discussed,  however. 

The  STA6SC-1  method  of  plasticity  calculations  between  MNR  iterations  appears 
therefore  to  raise  some  questions  about  its  use  in  dynamic  plastic  or  static 
cyclic  loading  problems. 

The  STAGSC-1  implementation  of  MNR  seems  to  be  basically  efficient  in  that  the 
direct  calculation  of  the  vector  {L ( x ) }  (eqs.  4.4,  4.5)  is  analogous  to 
the  method  of  calculating  a  pseudo-force  vector  to  account  for 
nonlinearities.  Moreover,  the  strategy  parameters  available  to  the  user 
provide  a  degree  of  control  over  the  MNR  procedure  which  is  not  available  in 
other  programs,  e.g.,  AD1NA  [11,12]. 

4.6  SOLUTION  OF  EQUATIONS 

STAGSC-1  employs  a  conventional  Cholesky  triangular  decomposition  with  forward 
and  backward  substitution  for  solution  of  equations.  Storage  is  based  on  the 
"skyline"  vector  concept  in  which  no  zero  elements  beyond  the  last  non-zero 
element  in  a  row  are  stored.  The  skyline  vector  stores  the  location  of  the 
last  non-zero  element  in  a  given  row.  The  procedure  is  outlined  in  the 
theoretical  manual  [3]  and  discussions  may  be  found  in  the  literature,  e.g.. 
Bathe  and  Wilson  [13].  No  other  solution  options  are  currently  available  in 
the  program. 
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4.7  FIGENVALUE  ANALYSIS 


Eigenvalue  analysis  is  required  for  both  buckling  and  vibration  analyses.  In 
principle,  there  is  no  fundamental  difference  between  the  two  eigenvalue 
problems.  However,  in  practice,  the  bifurcation  buckling  problem  requires 
solution  normally  only  for  the  smallest  buckling  load.  For  certain  types  of 
buckling  (e.g.,  pure  membrane  shear)  there  can  be  eigenvalues  which  are  equal 
but  opposite  in  sign.  For  vibration  problems,  the  eigenvalues  must  always  be 
positive  and  a  large  set  of  eigenvalues  and  eigenvectors  may  need  to  be 
determined.  Thus,  there  is  a  need  to  employ  a  method  which  is  suitable  for 
both  types  of  problem  and  which  is  relatively  "rugged,"  i.e.,  capable  of 
yielding  satisfactory  solutions  for  a  wide  variety  of  modeling  situations. 

The  method  implemented  in  STAGSC-1  is  basically  the  subspace  iteration 
method.  This  is  described  in  depth  by  Bathe  and  Wilson  [13]  and  somewhat 
sketchily  in  the  theoretical  manual  T3].  Curiously,  thp  manual  does  not  state 
explicitly  that  this  is  the  technique  being  used. 

Subspace  iteration,  as  described  in  Ref.  [13]  simultaneously  obtains  a  reduced 
number  of  eigenvectors.  An  initial  choice  is  made  of  a  set  of  preliminary 
independent  vectors  X  which  are  said  to  span  a  subspace  of  the  complete  set  of 
M  eigenvectors.  A  single  inverse  iteration  step  is  performed  in  which  a  new 
set  X  is  obtained  from  solution  of  the  equation 

K7  =  MX,  4.6 

where  both  K  and  M  are  of  order  mxm  and  X  and  J  are  of  order  mxp.  A  new 
eigenvalue  problem  is  then  solved  in  terms  of  the  reduced  matrices  K  and  W 


which  are  obtained 

from 

, 

and  K=rKX 

(pxp) 

4.7 

R  =  x 1 1  M  X  . 

(pxp) 

4.8 

The  reduced  eigenvalue  problem  can  be  stated  as  KQ  =  A  MQ 

4.9 

where  A  is  the  matrix  of  eigenvalues. 
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An  improved  estimate  of  the  subspace  vectors  X  is  then  obtained  as 


X  =  XQ  4.10 

at  which  point  the  whole  process  is  then  repeated. 

STAGSC-1  implements  this  method  with  some  modifications.  As  is  pointed  out  by 
Bathe  [13],  several  inverse  iteration  steps  may  be  performed  before  solving 
the  reduced  eigenvalue  problem  (in  order  to  reduce  the  number  of 
eigensolutions  performed).  This  is  done  in  STAGSC-1,  with  orthonormalization 
of  the  subspace  vectors  with  respect  to  K  or  M  at  each  step  of  the  iteration. 
Normally,  two  inverse  iterations  are  performed  at  the  start  and  three  (per 
eigensolution)  after  the  first  reduced  eigensolution.  Two  other  extra 
features  are  incorporated.  The  first  is  the  introduction  of  a  spectral  shift 
parameter  o  in  the  inverse  iteration  sweep.  Thus  the  solution  obtained  is 
for  the  equations 


[K  -  0M]  X  =  MX  4.11 

and  the  eigenvalues  obtained  for  the  reduced  system  are 

A*  =  A  -  ol  4.12 

The  second  modification  is  to  accelerate  the  convergence  of  the 
orthonormal i zed  vectors  by  means  of  Chebyshev  polynomials.  This  step  is 
performed  before  the  solution  of  the  reduced  eigensystem.  According  to  the 
developers,  this  technique  has  been  known  to  be  not  always  effective, 
particularly  with  some  computer  installations.  Its  use  is  controlled 
internally  and  depends  on  the  convergence.  Figure  4.3  provides  a  qualitative 
flow  chart  of  the  major  functions  in  the  eigensolution  system.  The  method  of 
solving  the  reduced  eigenvalue  problem  is  a  combination  of  Householder's 
tridiagonal ization  transformations  with  the  QR  method  of  extracting 
eigenvalues. 
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4.8  TRANSIENT  INTEGRATION 


Transient  dynamic  analysis  for  linear  and  nonlinear  shell  structures  is  one  of 
the  major  capabilities  in  STAGSC-1.  The  numerical  integration  with  respect  to 
time  of  the  equations  of  motion  can  be  accomplished  for  general  time  histories 
of  applied  load  or  displacement.  The  program  contains  a  library  of  five 
different  algorithms  for  performing  the  transient  integration  which  consists 
of  one  explicit  scheme  (central  difference)  and  four  implicit  methods.  The 
implicit  schemes  are  as  follows: 

A.  Trapezoidal  (Newmark,  p  =  1/4) 

B.  Gear's  second  order 

C.  Gear's  third  order 

D.  K.C.  Park's  method 

It  is  clear  at  the  outset,  that  for  such  an  array  of  options  for  implicit  time 
integration  to  be  useful,  the  program  user  needs  to  possess  a  greater  than 
average  level  of  sophistication  in  order  to  make  an  appropriate  choice.  The 
users  manual  [2]  seeks  to  minimize  this  difficulty  by  generally  recommending 
the  trapezoidal  rule  for  linear  structures  and  the  Park  method  for  nonlinear 
structures.  The  theoretical  manual  [3]  discusses  the  background  of  the  central 
difference  (explicit)  method  as  well  as  the  implicit  methods.  In  order  to  aid 
in  the  discussion  of  these  methods.  Figure  4.4  has  been  provided  in  an  attempt 
to  make  clear  the  basic  differences  between  the  methods  available  in  STAGSC-1 
and  other  methods  in  common  use.  The  characteristics  of  the  STAGSC-1  methods 
with  regard  to  stability,  numerical  damping  and  frequency  distortion 
(dispersion)  are  summarized  in  Table  4.3.  Although  the  central  difference 
method  has  been  implemented,  it  is  not  particularly  well  suited  to  its  use  in 
shell  problems  where  lower  frequency  modes  usually  dominate  the  response  of 
the  system.  The  low  stability  limit  requires  a  very  small  time  step  which  can 
largely  offset  the  inherent  efficiency  of  the  explicit  approach.  For  the 
majority  of  problems  therefore,  the  concern  must  be  with  the  relative  merits 
of  the  four  implicit  methods  provided.  With  the  exception  of  Gear’s  3rd  order 
method  (G-j),  the  implicit  methods  are  unconditionally  stable  for  linear 
conservative  systems.  This  is  referred  to  in  the  literature  as  A-stability 
[15,16].  G^  has  a  very  small  region  of  instability  for  a  combination  of 
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small  time  steps,  low  frequency  and  low  damping.  For  non-conservative 
systems,  Gpt  G3  and  the  K.C.  Park  method  (KCP)  are  conditionally  stable  as 
are  the  Wilson-e  and  Houbolt  methods.  The  trapezoidal  method  is  not  stable 
under  these  circumstances. 

The  question  of  greater  interest  to  a  potential  user  of  STAGSC-1  is  the 
performance  of  these  methods  in  the  context  of  nonlinear  response.  The 
stability  of  time-integration  operators  for  nonlinear  structural  dynamics 
problems  has  also  been  discussed  in  the  literature  (e.q..  Refs.  If  and  17). 
Implicit  operators  which  are  unconditionally  stable  for  linear  problems  have 
been  observed  to  exhibit  instability  in  some  nonlinear  problems.  The 
inference  has  been  drawn  that  the  stability  properties  of  the  implicit 
operators  are  lost  or  modified  in  the  nonlinear  regime.  Reference  17 
(Belytschko  an  Schoeberle)  presents  an  energy-based  proof  of  unconditional 
stability  for  the  Newmark-a  method  (b=1/4)  for  the  case  of  material 
nonlinearity.  The  proof  is  subject  to  the  restriction  that  the  internal 
energy  must  increase  monotonical ly  with  strain  and  remain  positive-definite. 
The  assumption  is  made  that  the  unconditional  stability  is  preserved  when 
geometric  nonlinearities  are  also  present  provided  that  the  requirements  on 
the  internal  energy  are  still  met.  It  is  concluded  that  the  loss  of  stability 
in  some  applications  is  due  to  errors  accumulated  during  the  solution  process 
and  not  to  the  integration  operator  per  se. 

A  somewhat  different  conclusion  is  arrived  at  in  Reference  15  (Park)  in  which 
nonlinear  stability  equations  for  a  number  of  implicit  operators  are 
developed.  These  criteria  apparently  include  the  character  of  the 
nonlinearity  (e.g.,  hardening  or  softening).  Stability  boundaries  are 
obtained  for  the  Houbolt,  Wilson,  Park  and  Newmark-a  methods.  It  is 
concluded  that  the  approximations  inherent  in  the  solution  procedures 
(initial-strain  or  tanqent  stiffness  methods)  are  responsible  for  the 
departures  from  unconditional  stability. 

Thus,  in  a  sense,  both  evaluations  arrive  at  the  same  conclusion  ( i . e . , 
unconditional  stability  is  affected  by  solution  approximations)  but  the 
implications  are  quite  different. 
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It  would  appear  that  on  the  basis  of  stability  the  methods  provided  are 
probably  quite  satisfactory  for  nonlinear  analysis,  but  the  user  ought  to  be 
aware  that  any  given  problem  should  be  evaluated  with  respect  to  its  likely 
nonlinear  behavior  when  choosing  the  method  of  integration.  With  respect  to 
accuracy,  the  four  implicit  methods  are  comparable.  It  does  appear,  for  a 
nonlinear  response  that  the  damping  and  dispersion  of  KCP  are  better  than  the 
Houbolt  method. 

4.9  USER  CODED  SUBROUTINES 

As  may  be  expected  in  a  program  developed  for  the  solution  of  nonlinear 
problems,  STAGSC-1  has  a  significant  capability  for  the  user  to  provide  his 
own  codinq  for  problem  specification  and  execution.  Thus,  there  are  a  total 
of  thirteen  dummy  subroutines  in  the  program  for  which  the  user  can  provide 
FORTRAN  coding.  Some  of  the  subroutines  provide  additional  capability  while 
others  are  mainly  used  to  reduce  the  bulk  of  input  data.  Each  of  the 
subroutines  will  be  briefly  described  and  commented  on. 

A.  CROSS--This  routine  defines  cross-section  dimensions  and  material 
properties  for  beams  and  stiffeners.  Geometry  and  material  type  can 
be  specified  as  functions  of  spatial  coordinates.  Used  in  addition 
to,  or  instead  of  data  cards. 

B.  FORCET--Descr ibes  variation  of  load  factor  with  time  for  load  system 
A  or  B.  Used  instead  of  data  cards. 

C.  UGRID— Allows  independent  specification  of  grid  coordinates  for  mesh 
generation  in  terms  of  reference  surface  geometry.  Necessary  for 
quadrilateral  elements  in  a  quadrilateral  plate. 

D.  LAME — Allows  definition  of  a  shell  unit  geometry  not  included  in  the 
twelve  built-in  options. 

E.  SKEWS — Defines  orientation  on  shell  surface  of  the  attachment  line  of 
a  discrete  stiffener  which  does  not  follow  the  reference  surface 

gr idl ines. 
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F.  TEMP--This  is  the  only  means  by  which  temperatures  can  be  specified. 
Temperatures  can  vary  with  surface  coordinates  and  through  a 
stiffener  cross-section. 

G.  UCONST--Def ines  linear  constraint  conditions  using  Lagrange 
multipliers.  Must  not  be  used  with  explicit  integration. 

H.  UPRESS — Defines  space  and  time  variation  of  pressure  loads.  Pressure 
may  be  “follower"  or  "live"  loading. 

I.  USRLD--Def ines  spatial  variation  of  loads  including  initial 
displacements  and  velocities.  Saves  preparation  of  bulky  sets  of 
data  cards. 

J.  USRELT--Def ines  element  connectivity  for  an  element  unit.  Essential 
when  there  are  more  than  a  few  elements. 

K.  USRPT--Def ines  node  point  geometry  for  element  units.  Essential  when 
there  are  more  than  a  few  nodes.  May  also  define  additional  points 
in  a  shell  unit. 

L.  WALL--Def ines  shell  wall  construction  (layers,  composite  material, 
stiffeners,  wall  thickness)  and  material  properties  which  may  vary 
over  the  reference  surface.  Material  properties  may  not  vary  within 
a  layer.  Does  not  apply  to  plasticity  data. 

M.  WIMP--Def ines  small,  geometric  perturbations  of  the  shell  from  the 
reference  surface  in  terms  of  first  spatial  derivatives. 

4.10  RESTART  CAPABILITY 

For  practical,  nonlinear  structural  analysis,  a  useful  computer  program  must 
include  a  flexible  restart  capability.  Ideally,  the  user  should  be  able  to 
restart  the  analysis  at  any  desired  point  so  that  a  different  loading  strategy 
may  be  used  or  perhaps  an  eigenvalue  solution  obtained.  STAGSC-1  has  such  a 
capability.  At  the  user's  option,  three  separate  files  may  be  saved  for 
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restart.  (TAPE22,  2 3  and  24).  The  solution  is  saved  on  TAPF22  while  TAPE  23 
and  TAPE24  contains  the  stiffness  matrix  and  the  factored  stiffness  matrix. 
Normally,  upon  restart,  refactorization  occurs  but  the  user  can  override  this 
by  using  TAPE24.  The  contents  of  TAPE22  include  displacements,  velocities  and 
plastic  strains,  depending  on  the  type  of  analysis.  Additionally,  stresses, 
strains  and  stress  resultants  can  be  saved  on  the  same  file  for 
post-processing. 

STAGSC-1  provides  the  option  to  save  data  either  at  every  load  (time)  step  or 
from  the  final  three  steps.  This  could  be  improved  by  permitting  saving  at 
specified  load  step  intervals  as  in  the  MARC  program.  The  advantage  of  this 
is  that  the  flexibility  to  restart  at  a  number  of  stages  is  retained  but  with 
substantial  savings  in  file  space.  This  can  be  a  significant  factor  for  the 
analysis  of  a  real  nonlinear  problem. 

4.11  INPUT  AND  OUTPUT 

Input  and  output  are  often  the  basis  for  user  attitudes  towards  a  structural 
analysis  program.  Factors  which  influence  these  attitudes  are  many  but  the 
major  ones  are 

o  logical  input  flow 

o  input  format  (free  form  or  otherwise) 
o  ability  to  provide  comments 

o  ease  of  generating  repetitive  data 

o  understandabi 1 ity  of  input  instructions 
o  control  over  output 

o  format  and  labeling  of  output 

A  user's  reaction  to  the  input  required  for  a  new  program  is  often  influenced 
by  experience  with  other  programs  which,  of  course,  may  place  the  new  program 
in  a  good  or  bad  light  depending  on  the  previous  experience.  Nevertheless, 
the  concept  has  arisen  of  "user  friendliness"  as  a  measure  of  the  attitude 
which  a  program  may  develop  in  the  user.  The  meaning  of  this  is  obviously 
subjective,  but  in  such  a  context  STAGSC-1  would  probably  rate  as  average. 

Ease  of  input  depends  strongly  on  the  logical  flow  of  the  input  stream  and  the 
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understandabil ity  of  the  input  instructions.  STAGSC-1  input  is  quite  logical 
and  if  the  input  manual  is  followed  carefully  a  new  user  may  expect  to  obtain 
a  set  of  input  which  will  execute  successfully  after,  perhaps,  a  couple  of 
tries.  However,  it  is  apparent  that  both  the  input  logic  and  the  input 
instructions  are  strongly  linked  to  the  programming  and  are  not  based  on  some 
concept  of  what  might  constitute  "good"  input.  This  is  not  necessarily  a 
criticism  of  the  input  but  more  a  description  of  the  type  of  input.  Excellent 
features  are  the  free-form  input  and  the  ability  to  include  user  commenting. 
STAGSC-1  input  would  lend  itself  eadily  to  the  interactive  mode  since  the 
instructions  used  in  the  manual  are  already  in  the  required  form. 

If  the  user  selects  any  of  the  twelve  standard  shell  units,  the  input  required 
for  the  generation  of  bulk  data  is  minimal.  Also,  the  constraints  that 
provide  compatibility  between  shell  units  are  easily  imposed.  For  other 
geometries,  the  user  written  subroutine  LAME  may  be  used  to  define  the 
reference  surface.  For  element  units,  the  choice  is  either  individual  input 
of  each  node  and  element  or  automatic  generation  using  user  subroutines  USRPT 
and  USRELT.  Individual  node  and  element  input  is  unacceptable  for  more  than  a 
few  elements  so  the  use  of  subroutines  is  almost  always  required.  An 
improvement  would  be  the  inclusion  of  simple  linear  mesh  generators  and 
element  pattern  generators  to  provide  a  rapid  means  of  generating  meshes  for  a 
large  class  of  problems. 

For  load  input,  the  ability  to  input  both  concentrated  and  distributed  loads 
either  individually  at  mesh  points  or  along  specified  rows  and  columns  is  a 
good  basic  feature.  This,  toqether  with  the  user  subroutines  USRLD  and  IIPRESS 
for  generating  loads  provide  a  generally  satisfactory  capability.  For  element 
units,  not  all  load  input  options  are  yet  operational,  e.g.,  distributed 
forces  and  moments. 

Output  control  may  be  exercised  separately  on  displacements,  strains, 
stresses,  stress  resultants,  stresses  and  strains  at  yielded  points  and 
forces.  The  data  for  each  shell  unit  is  output  as  a  block.  However,  the 
frequency  may  be  specified  differently  for  each  output  quantity 
(displacements,  stresses,  etc.)  and  also  for  each  shell  unit.  Based  on  some 
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limited  experience,  it  seems  that  if  different  frequencies  are  specified  for 
say  stresses  and  displacements,  both  will  be  printed  out  at  the  higher 
frequency  which  is  contrary  to  the  manual. 

In  addition,  selected  Stresses  and  displacements  may  be  output  at  each  load  or 
time  step.  Thus,  if  it  is  desired  to  output  data  for  a  certain  portion  of  the 
structure  at  certain  load  step  intervals,  the  only  means  of  doing  this  is  to 
specify  a  special  shell  unit  just  for  this  region  of  the  structure  (since  all 
displacements,  etc.,  are  output  for  a  shell  unit).  This  may  be  inconvenient 
and  therefore  some  additional  selectivity  is  needed  to  cater  for  such  a 
situation. 

4.12  POST-PROCESSING  AND  PLOTTING 

Currently,  post-processing  with  STAGSC-1  is  only  partially  operational.  Its 
developmental  status  is  not  yet  comparable  with  that  of  the  analysis  program 
which  is  quite  a  serious  disadvantage  when  performing  nonlinear  analysis. 
Moreover,  in  the  case  of  a  nonlinear  analysis,  the  user  often  needs  the 
ability  to  access  the  solution  data  file  and  perform  his  own  post-processing 
directly.  A  typical  requirement  would  be  to  extract  inelastic  strain 
histories  from  the  solution  and  post-process  them  according  to  design  code 
criteria,  e.g.,  the  ASME  Boiler  and  Pressure  Vessel  Code,  Code  Case  N-47.  The 
necessary  descriptions  of  the  structure  of  the  solution  files  are  not 
available  in  the  documentation  however,  so  this  option  is  unavailable  to  the 
user. 

The  STAGSC-1  post-processing  program  (STAPL)  is  executable  in  tandem  with  the 
analysis  or  separately  by  saving  the  post-processing  file  (TAPE22).  STAPL  has 
been  developed  from  routines  published  by  NASA  [18]  which  provide  deformed  and 
undeformed  geometry  plots  and  also  contour  plots.  The  separate  routines  were 
merged  by  Anamet  Laboratories  [19]  and  further  developed  by  Lockheed.  The 
current  range  of  capabilites  listed  in  the  user's  manual  [2]  are  as  follows: 

A.  geometry  plots--deformed,  undeformed  and  exploded  (useful  for  mesh 
checking) 
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B.  contour  plots--displacements,  velocities,  force  residuals*, 

eigenvectors,  force  resultants,  strains,  stresses,  temperatures*, 
masses*,  initial  conditions  and  loads. 

The  items  marked  with  an  asterisk  are  not  yet  operational  accordinq  to  the 
manual.  Other  features  advertised  which  are  also  not  operational  are  geometry 
plots  showing  only  shell  unit  boundaries,  solution  quantities  as  vectors 
emanating  from  the  nodes;  and  automatic  plot  scaling,  (valuation  of  STAGSC-1 
plottinq  is  not  included  in  the  scope  of  this  report. 

Other  development  work  has  included  the  interfacing  of  the  GIFTS*  interactive 
graphics  package  for  mesh  generation  with  STA6SC-1.  This  combination  is 
currently  being  evaluated  by  ONR. 

4.13  SPECIAL  STAGS  FEATURES 

This  section  is  intended  to  highlight  those  features  of  the  STAGSC-1  program 
which  serve  to  differentiate  it  from  other  finite  element  nonlinear  structural 
analysis  programs.  To  begin  with,  STAGSC-1  is  the  outcome  of  approximately 
thirteen  years  of  development  effort  in  an  aerospace  environment.  It  is  this 
environment  which  has  stimulated  the  development  of  a  number  of  analysis 
programs  for  shells  of  revolution,  e.g.,  BOSOR,  DYNAPLAS,  SATANS,  etc.  [20]. 
STAGSC-1  is  an  outgrowth  of  this  effort  which  extends  the  capabilities  to 
general,  three-dimensional  thin  shells.  Other  finite  element  analysis 
programs  (e.g.,  MARC)  have  incorporated  shell  elements  in  their  element 
libraries  but  supposedly  cannot  match  the  greater  efficiency  of  the  special 
purpose  shell  programs. 

Therefore,  the  basic  advantage  of  STAGSC-1  is  its  emphasis  on  shell  analysis. 
The  inclusion  of  beam  and  spar  type  elements  does  not  make  it  a  natural  choice 
for  solely  beam  or  truss  types  of  structure,  although  it  certainly  is  able  to 


♦GIFTS  is  a  finite  element  mesh  generation  and  analysis  program  originally 
developed  for  the  analysis  of  ship  structures  and  supported  by  ONR. 


0918B-84B.-2 
( S3034)  24 


59 


perform  such  analyses.  In  this  context  one  may  pick  out  a  number  of  features 
of  the  program  which  are  probably  unique  to  STA6SC-1.  These  will  be  discussed 
briefly  in  the  remainder  of  this  section. 


A.  Shell  Unit  Concept 

This  has  already  been  introduced  and  described  in  Section  4.2.  The 
distinguishing  feature  is  the  underlying  grid  on  which  the  element 
mesh  can  be  overlaid.  This  makes  it  possible  to  define  a  library  of 
standard  geometries  with  a  minimum  of  input.  Thus  the  mesh 
variability  which  may  be  required  for  a  specific  problem  can  be 
divorced  from  the  generation  of  the  surface  geometry. 

B.  Initial  Imperfections 

The  ability  to  specify  an  imperfect  geometry  as  small  perturbations 
to  a  basic  reference  geometry  is  a  feature  of  great  value  in  a 
program  oriented  towards  shell  buckling  and  collapse.  To  achieve 
this  in  a  general  purpose  program,  if  possible  at  all,  would  probably 
require  the  writing  of  a  special  mesh  generator  for  each  problem.  In 
STAGSC-1 ,  even  if  the  reference  surface  is  not  part  of  the  library, 
the  writing  of  the  LAME  subroutine  for  the  reference  surface  and  WIMP 
for  the  initial  imperfections  is  likely  to  be  the  most  convenient  way 
of  generating  the  data. 

C.  Layered  and  Composite  Shell  Wall  Construction 

Relatively  complicated  shell  wall  designs  can  be  handled  by 
STAGSC-1.  The  types  which  may  be  included  are  as  follows: 

(1)  multiple  anisotropic  layers 

(2)  multiple  fiberwound  layers 

(3)  walls  reinforced  by  corrugated  skin 

(4)  wall  properties  defined  by  a  shell  wall  stiffness  matrix 

(5)  any  of  the  above  with  "smeared"  stiffeners 

b() 
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Multiple  anisotropic  layers  are  specified  by  layer  thickness, 
orthotropic  elastic  material  properties  and  a  principal  direction  for 
these  properties.  This  can  also  be  used  for  fiberwound  materials 
where  the  layer  properties  are  available.  Corrugated  skin 
reinforcement  is  modeled  by  trapezoidal  shaped  corrugations  whose 
cross-section  dimensions  are  input.  The  shell  wall  stiffness  matrix 
method  is  available  for  layered  composite  walls  whose  overall 
stiffness  properties  are  known. 

0.  White-Bessel ing  Plasticity  Model 

This  is  not  a  commonly  implemented  constitutive  model  of  plasticity 
although  the  Mroz  model,  to  which  it  is  related,  is  available  in  the 
PLANS  program  [20].  The  user  should  be  aware  that  although  there  is 
some  evidence  [21]  to  suggest  that  methods  based  on  the  mechanical 
sublayer  concept  model  reversed  loading  behavior  well  for  some 
materials,  there  is  not  as  yet  a  substantial  body  of  testing  or 
analytical  experience  to  validate  its  use  fully.  The  potential  user 
should  therefore  be  prepared  to  perform  his  own  validation  for  his 
application. 

E.  Library  of  Load-Time  Histories 

STAGSC-1  contains  three  specific  load-time  histories  for  transient 
integration.  These  are 

(1)  trapezoidal  variation 

(2)  trignometric  variation--sinusoidal ,  cosine  square  impulse  or 
cosine  square  ramp  functions 

(3)  linear  ramp  and  exponential  decay 

These  basic  ingredients  can  model  a  significant  variety  of 
transients.  However,  more  general  forcing  functions  may  require  the 
user-supplied  subroutine  FORCET. 


Cl 
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TABLE  4.1 

SHELL  SURFACE  GEOMETRIES 


No.  of 


SHELL 

Description 

Coordinates 

Comments 

2 

Rectangular  Plate 

4 

4  edge  coordinates 

3 

Quadrilateral 

Plate 

8 

8  corner  coordinates 

4 

Annular  Plate 

4 

2  radii,  2  subtended  angles 

5 

Cyl inder 

5 

2  axial,  2  subtended  angles,  1 
radius 

6 

Cone 

6 

2  axial,  2  subtended  angles,  2 
radii 

7 

Sphere 

5 

2  meridional  and  2  azimuthal 
angles,  1  radius 

8 

Torus 

6 

2  meridional  and  2  axial  angles, 
bend  radius,  cross-section  radius 

9 

Ell iptic  cone  or 
or  Cylinder 

7 

2  axial,  2  angles,  2  major  and 

1  minor  radii  (similar  cross- 
sections) 

10 

Paraboloid 

6 

2  axial,  2  angles,  distances  from 
apex  to  focus  and  to  smaller  end 

11 

Ellipsoid 

6 

2  meridional  and  2  azimuthal 
angles,  major  and  minor  radii 

12 

Hyperboloid 

7 

2  axial,  2  angles,  3  coordinates 
to  define  asymptote 
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TA8LE  4.2 

SUMMARY  OF  ELEMENT  LIBRARY 


I 

I 

I 


cn 

c 

1 

1 

cj 

o 

D 

u 

i 

1 

1 

1 

1 

■r- 

•F“ 

•r* 

•1“ 

i 

1 

1 

to 

3 

1 

1 

XI 

JO 

jO 

3 

i 

1 

1 

c 

C 

t 

1 

3 

3 

3 

3 

i 

1 

1 

o 

C V 

1 

1 

<_> 

O 

O 

C_> 

» 

1 

1 

T- 

CO 

4-1 

u 

c 

3 

LX. 

C 

c 

CU 

<D 

o 

o 

31 

q. 

•r— 

•r— 

3 

rra 

4-> 

c 

CU  <U 

_c 

rd 

rd 

o 

CD 

L/> 

r— 

< — 

•r— 

O  3 

<J\ 

on 

-M  -»-> 

4->  CU 

l+_ 

C 

C 

00  rd 

o 

rd 

•d 

•r»  t- 

O 

■K  O 

cu 

4->  S- 

4-»  S- 

i  00 

«r- 

—  +■> 

c 

t/>  4-> 

00  +-> 

+J  c 

+-> 

— 

<u 

n 3 

S- 

*r~ 

s. 

rd 

rd 

+ 

3 

rd 

5  u 

rd 

S  o 

3  S- 

s- 

CJ  s- 

o 

u 

Q_ 

<D 

4->  *»- 

cu 

■M  *r- 

•r-  4-> 

3 

♦r- 

•r- 

o 

I 

C 

X-> 

c 

4-> 

4-> 

rd 

4-J 

4-> 

c 

•r* 

L_  rd 

•r— 

i-  rd 

<d  3 

3 

rd 

rd 

t— ( 

f— 

rd  i- 

r— 

c0  S- 

l-  rd 

cr 

rd 

S-  rd 

i- 

QJ  3 

a; 

3 

a> 

3  QJ 

3 

, — 

c  rd 

, — 

c  <d 

rd  C 

t — 

c 

rd  C 

rd 

•r—  3 

r— 

•r-  3 

3  *»- 

r— 

•r* 

3  *i— 

3 

<c 

—1  cr 

<c 

-3  O’ 

O'  -J 

-3 

O’  _j 

O' 

x:  • 

4-» 

• 

u  i- 

O 

»d 

on 

on 

rd  CD 

rd 

*o 

3 

(U  4-» 

a ) 

3 

to 

c 

c  • 

c 

• 

C 

oo  <U 

to  4~> 

a; 

a;  v- 

4-> 

■m  a> 

4->  S- 

4-J  d 

3  3 

3  C 

<v 

*d 

<d  U 

rd  CU 

rd 

r—  *r— 

i —  cu 

4-> 

4->  4-> 

■M 

c 

Q.  00 

CL  E 

rd  c 

00 

cn  4-> 

VI  c 

oo  O 

• 

3 

CU 

cu 

c 

C  rd 

c  cu 

C 

£-  *r- 

S-  u 

4-> 

4_>  Q 

o 

O 

o  u 

O  +-» 

cu 

0J  E 

CU  rd 

tO 

</> 

•f— 

*•—  4-» 

•r-  <d 

c 

C 

C  ' — 

•r— 

■r-  +-> 

4-> 

■M  C 

-M  4-> 

4->  •— 

s- 

-M 

S-  CL 

3:  rd 

<d 

rd  QJ 

ro  rd 

rd  00 

o 

O  <d 

o  CO 

4-> 

-4-1 

4-» 

4->  E 

4-> 

4-1  C 

u 

U 

CJ  •»“ 

• 

x-> 

o 

O  OJ 

o  c 

o  rd 

4-> 

3 

t 1 

-o 

3  c 

S- 

u 

V-  O 

^  s- 

JZ 

^  C 

JZ 

c 

c  cu 

rd 

*r— 

4-> 

u 

U  O) 

u  *— 

o 

rd 

rd  £ 

co 

CO  1— 

ro  -m 

CO 

rd 

*P  E 

rd  rd 

0/ 

Q. 

rd 

r-  • 

CU 

CU  CU 

(U  *r- 

o 

4-> 

4~)  U 

■a 

*3  l/l 

3  -*-> 

3  <d  k. 

u 

4-» 

c 

c  ro 

c 

C  -r- 

C  O 

C  •»-  CU 

4-) 

4-1  rd 

4-1  C 

(4_ 

a> 

a)  r— 

<d 

rd  *3 

rd 

rd  X  4-> 

rd 

rd  r— 

rd  CU 

o 

£ 

E  Q- 

rd  C 

Q. 

31 

s 

cu  to 

to 

00  i — 

oo  i — 

on  CU 

to 

to  oo 

tO  3  • 

OJ 

<j 

<j  •<— 

c 

C  rd 

C  rd 

C  CO  O 

c 

C  -r- 

C  rd  OO 

Q. 

rd 

*d  *0 

o 

o  •*- 

O  -r- 

O  3 

o 

O  3 

0  4-10) 

>», 

«r— 

*»“  X 

•r-  X 

•r-  r—  +-» 

•r- 

•f 

*»—  3 

CL 

r>  r— * 

4-> 

■fJ  rd 

4->  rd 

3  Q.  <d 

4-> 

4->  •— 

4-13  0 

(/I 

i/>  rd 

<d  • 

rd 

rd 

rd 

rd 

rd  rd 

rd  C  C 

•r" 

«—  3 

• —  00 

i —  OO 

r-  3  C 

r-“ 

r— “  *f“ 

r—  rd 

•o 

3  X 

00  c 

00  3 

to  3 

V>  c  o 

to 

0O 

to  0) 

rd 

C  CL) 

C  r— 

C  r— 

C  <U  *»- 

c 

C  C  • 

C  r—  3 

rd 

rd  Cl 

rd  CL 

rd  4-> 

rd 

rd  CU  to 

rd  rd  *r- 

to 

03  1/1 

i.  jt 

S- 

$~ 

l  r  re 

i- 

V-  31  <U 

w  E  co 

•r* 

•r-  3 

x->  o 

-«->  "O 

4-»  3 

u  uu 

4-> 

+■>  C  3 

4-1  t.  3 

>c 

X  ’ — 

<d 

c 

C 

rd  O 

rd  O 

o  •*— 

c* 

<c  o. 

ro  <u 

CO  <D 

oo  cu 

ro  cu  i- 

CNJ 

C\J  4->  C 

csj  c  E 

4-  • 

O  U_ 

•  O 

o  • 

z  o 


cu 

Q- 

>> 


c 

cu 

e 


V- 
« 3 
CO 


O 

o 

CNJ 


tr> 

CNJ 

CO 

CO 

VO 

Oi 

CM 

i 

i 

i 

i 

i 

1 

1 

1 

1 

f 

1 

1 

1 

i 

i 

i 

i 

i- 

rd 

*—  0) 

3  C 

31  rd 

i 

i 

i 

i 

f 

I 

X 

1 

cu 

z 

C  l- 

z 

s 

”i 

i 

CO 

t 

i 

1 

1 

I 

l 

rd  -O 

E 

l 

l 

1 

1 

» 

i 

i 

i 

1 

1 

t 

I 

^  a> 
h-  s: 

i 

l 

1 

1 

o 

r— * 

O 

r— 

o 

r— 

CJ 

o 

r— 

CJ 

CV) 

o 

o 

o 

CJ 

CJ 

CJ 

CJ 

CJ 

<n 

OC 

cn 

rd 

O. 


63 


Lr  =  perpendicular 
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Figure  4. 1(a)  Linear  Pre-buckling 


Figure  4. 1(b)  Nonlinear  Pre-buckling  Figure  4. 1(c)  Nonlinear  Pre-buckling 

-  Softening  -  Stiffening 


Figure  4.2  Nonlinear  Solution  Procedure 
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Figure  4.3  Subspace  Iteration  Algorithm 
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5.0  VERIFICATION  EXERCISES 


The  verification  procedure  followed  for  the  purposes  of  this  evaluation  was 
simply  to  execute  three  different  problems  which  were  part  of  the  STAGSC-1 
program  file  supplied  by  Lockheed.  The  purpose  of  the  verification  was  solely 
to  ensure  that  the  version,  as  modified  for  the  CDC  7600,  was,  in  fact, 
functioning  correctly  on  the  Westinghouse  system.  This  does  not  constitute 
verification  in  the  normal  sense*,  but  was  considered  to  be  sufficient  in  view 
of  the  fact  that  the  problems  to  be  run  in  the  course  of  the  evaluation  would 
themselves  provide  substantial  verification.  The  verification  problems  and 
their  solutions  will  be  described  in  this  section. 

5.1  SHALLOW  ARCH  PROBLEM 

The  first  verification  exercise  was  a  static,  nonlinear  analysis  of  a  shallow 
arch  with  a  uniform,  radially  inward  pressure  loading.  Figure  5.1  shows  the 
model  data.  The  boundary  conditions  are  simple  supports  at  each  end  of  the 
arch.  The  solution  converged  using  three  load  increments  with  6  to  8 
iterations  required  for  each  increment.  One  refactoring  occurred  during  the 
final  increment.  Figure  5.2  shows  the  center  displacement  plotted  against 
pressure.  The  solution  obtained  agreed  identically  with  the  output  generated 
by  Lockheed.  No  analytically  based  solution  for  this  problem  was  available 
for  purposes  of  comparison. 

5.2  BUCKLING  OF  AN  ANISOTROPIC  FLAT  PLATE 

The  program  options  involved  in  this  problem  include  an  eigenvalue  analysis 
for  buckling,  with  material  axes  rotated  with  respect  to  the  global  axes  by 
45°.  The  bifurcation  buckling  load  is  124.6  lb/in.  Figure  5.3  shows  the 
plate  geometry,  the  applied  loading  and  the  buckling  mode.  No 
analytical ly-based  solution  was  available,  but  the  results  agreed  with  the 
Lockheed  output. 


♦Verification  is  nominally  the  responsibility  of  the  developer,  and  is  intended 
to  ensure  that  the  capabilities  provided  by  the  developer  are,  in  an  isolated 
sense,  functioning. 
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5.3  CIRCULAR  RING  WITH  INITIAL  VELOCITY 


Figure  5.4  shows  the  geometry  for  a  circular  ring  segment  with  a  uniform, 
radial ly-inward  initial  velocity.  Radially  symmetric  boundary  conditions  give 
rise  to  a  solution  which  predicts  purely  radial  motion  uniform  along  the 
segment.  Figure  5.5  shows  the  time  history  of  the  radial  deflection 
response.  The  measured  period  yields  a  frequency  of  952.4  Hz  which  compares 
with  the  exact  value  of  1000  Hz  [36]  (-4.8%  error).  The  results  agreed 
identically  with  the  Lockheed  output. 

These  three  checks,  when  first  attempted,  failed  to  execute  and  identified  a 
coding  error  which  was  corrected  after  consultation  with  Lockheed.  The 
subsequent  successful  executions  confirmed  that  the  program  was  performing 
correctly  on  the  Westinghouse  system. 
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Figure  5.1  Shallow  Arch  Problem 
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Figure  5.2  Shallow  Arch  Problem 
-  Results 


Buckling  of  an  Anisotropic  Plate 


Figure  5.4  Circular  Ring  With  Initial  Velocity 


e  5.5  Ring  with  Initial  Velocity  -  Results 


6.0  ADVANCED  EVALUATION  EXERCISES 


The  most  effective  and  convincing  form  of  evaluation  of  a  structural  analysis 
computer  program  is  the  obtaining  of  results  from  the  execution  of  special 
problems  designed  to  test  in-depth  the  various  options  of  the  program  both 
individually  and  in  combination.  This  basic  idea  was  developed  in 
considerable  detail  by  Nickel 7  [1].  It  was  made  clear  in  that  discussion  that 
these  are  not  conventional  verification  exercises  but  problems  chosen  and 
executed  in  order  to  study  specific  aspects  of  the  behavior  of  the  elements  in 
the  program  library  together  with  the  solution  algorithms.  The  rigorous 
evaluation  of  a  sophisticated  program  such  as  STAGSC-1  based  on  these 
guidelines  is  a  very  large  task.  The  scope  of  the  present  work  is  therefore 
limited  to  provide  an  in-depth  study  of  some,  but  not  all,  of  the  capabilities 
of  the  program.  The  areas  selected  for  the  advanced  evaluation  are  as  follows 

o  Element  convergence 
o  Eigenvalue  extraction 

o  Transient  integration 

o  Nonlinear  solution  algorithm 

These  will  be  described  and  the  results  presented  in  the  four  following 
subsections. 

6.1  ELEMENT  CONVERGENCE 

A  dual  approach  has  been  adopted  in  investigating  the  convergence  of  STAGSC-1 
elements.  First  of  all,  a  direct  convergence  study  has  been  made  by  grid 
refinement  of  two  mul t i -element  problems  which  have  well  documented 
solutions.  The  problems  selected  are: 

A.  a  cylindrical  roof  (or  barrel  vault)  simply  supported  at  each  end  and 
subjected  to  gravity  loading. 

B.  a  square  flat  plate  with  clamped  edges  and  uniform  pressure  (linear 
and  nonlinear) . 
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These  have  been  extensively  studied  [14,  22,  23]  using  both  analytical  and 
numerical  methods  of  solution.  Also  the  convergence  properties  of  other 
elements  for  the  cylindrical  roof  are  well  known.  The  second  approach  is  to 
determine  the  eigenvector  components  of  the  element  stiffness  matrix  and 
resolve  the  energy  content  for  simple  loading  cases  into  the  various 
eigenmodes.  The  spectrum  of  eigerimodes  gives  insight  into  the  deformation 
responses  which  can  be  obtained  using  a  particular  element  and  can  also  show 
whether  there  are  any  spurious  (zero  energy)  modes  inherent  in  the  element 
formulation. 

In  this  work,  attention  has  been  confined  to  the  puadri 1 ateral  shell 
elements.  The  reasons  for  this  are: 

A.  STAGSC-1  is  primarily  a  shell  program; 
i  B.  The  triangular  shell  elements  are  embodied  in  the  quadrilaterals. 

i 

[  6.1.1  CYLINDRICAL  ROOF  PROBLEM 

The  qeometry  of  the  roof  is  shown  in  Figure  6.1.  For  gravity  loading,  there 
are  two  vertical  planes  of  symmetry  and  hence  only  one  quarter  of  the 
structure  needs  to  be  modeled.  The  majority  of  published  solutions  to  this 
problem  assume  a  zero  value  for  Poisson's  ratio. 

For  the  present  investigation,  the  two  classes  of  quadrilateral  plate  element 
(410  and  420  series)  were  tested.  Two  of  the  420  series  (420  and  421)  were 
used  but  both  failed  to  obtain  a  solution.  In  each  case  a  message  was  printed 
indicating  ill-conditioning  of  the  stiffness  matrix  and  termination  of  the 
solution  because  of  this.  Reasons  for  this  failure  are  not  understood  and  are 
beinq  studied  by  Lockheed.  The  mesh  configurations  used  for  the  1/4  model 
(see  Fiqure  6.1)  ranged  from  2  x  2  (4  elements)  to  10  x  10  (100  elements). 

Two  displacement  measures  of  convergence  were  used  and  gave  similar  results. 
These  were  the  vertical  deflections  at  the  mid-point  of  the  free  edge  (Point 
B)  and  the  mid-point  of  the  whole  roof  (Point  C).  Figures  6.2  and  6.3  show 
the  deflections  WD  and  W„  plotted  against  total  number  of  degrees  of 
freedom.  Fiqure  6.2  shows  also  the  theoretical  solution  for  deflection  Wg 
based  on  shallow  shell  and  deep  shell  theories  (see  Reference  23).  Probably 
the  best  numerical  solutions  to  this  problem  have  been  achieved  using  the 
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A8AQUS  [25]  and  the  MARC  [26]  finite  element  programs  and  these  are  included 
as  benchmark  comparisons.  In  addition,  the  shallow  curved  shell  triangular 
element  of  Cowper,  Lindberq  and  Olson  [24]  is  shown  as  representative  of  the 
earlier  approaches  to  shell  element  formulation  based  on  a  shallow  shell 
approx imat ion  (Novozhilov  theory).  The  ABAQUS  element  is  the  thin  shell 
i soparametr ic  element  due  to  Ahmad,  while  the  MARC  element.  24  is  based  on  the 
de  Veubeke  element  and  Koiter-Sanders  shell  theory  and  is  also  doubly-curved, 
isoparametric.  Both  of  these  elements  give  very  accurate  results  using  only 
four  elements. 

The  STAGSC-1  elements,  being  flat,  approximate  the  geometry  relatively  crudely 
for  the  coarsest  mesh  (4  elements)  and  the  results  are  correspondingly 
inaccurate.  The  410  element  for  the  4  element  mesh  overestimates  the  free 
edge  displacement  Wg  and  the  center  displacement  by  about  41%  in  each 
case.  The  411  element  performs  better  for  the  coarse  mesh  but  uses  99  degrees 
of  freedom  instead  of  63  (410).  However,  the  convergence  properties  of  the 
411  element  appear  rather  worse  than  the  410  despite  the  higher  order  shape 
functions.  Moreover,  its  initial  high  convergence  rate  gives  an  overshoot 
with  respect  to  the  exact  solution  which  is  why,  in  the  end,  it  does  not  give 
any  better  results  than  the  410  element.  The  NASTRAN  evaluation  study  [14] 
gives  a  very  comprehensive  comparison  between  the  NASTRAN  element  TRSHL  and  a 
number  of  other  elements  including  the  Ahmad  isoparametric.  Comparing  the 
STAGSC-1  results  with  this  data  it  is  obvious  that  the  STAGSC-1  elements 
perform  considerably  better  than  TRSHL  but  not  so  well  as  MSC/NASTRAN's  QUAD  4 
element . 

In  summary,  it  can  be  stated  that  the  410  series  performs  acceptably  for  a 
curved  shell  problem  bearing  in  mind  that  the  element  formulation  is  for  a 
flat  plate. 

6.1.2  FLAT  PLATE  PROBLEM 

The  failure  of  the  420  series  elements  to  provide  a  solution  to  the 
cylindrical  roof  problem  made  it  necessary  to  select  a  different  test  case.  A 
flat  square  plate,  simply  supported  or  clamped  at  the  edges  was  used  in  the 
NASTRAN  evaluation  study  by  Jones,  et.  al.,  [14]  to  obtain  convergence  data 
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for  the  QDPLT  (quadr i lateral  plate)  element.  A  purely  linear  solution  to  this 
problem  will  only  provide  data  for  the  bending  behavior  of  an  element. 

herefore,  in  order  to  test  both  bending  and  membrane  convergence,  the  large 
deflection  problem  must  be  solved.  Linear  solutions  are  available  in  Roark 
[27]  for  the  flat  plate  subject  to  a  wide  variety  of  loads  and  boundary 
conditions.  However,  for  the  nonlinear  case  the  only  solutions  provided  are 
for  uniform  pressure  loading  with  various  edqe  conditions.  The  case  selected 
therefore  was  the  square  plate,  clamped  along  all  four  edges  and  subject  to 
uniform  pressure  loadinq.  The  linear  solution  to  this  problem  in  terms  of  the 
deflection  at  the  mid-point  of  the  plate  is 


W.  =  0.01376  2^-,  where  A  denotes  the  plate  center. 

Ei 

Eiqure  6.4  shows  the  geometry  of  the  plate.  For  this  geometry  and  a  pressure 
of  12500  p.s.i.,  the  maximum  deflection  is  given  by 

=  0.2752  inches 

The  nonlinear  solution  is  presented  as  a  function  of  the  applied  pressure  and 
is  shown  in  Figure  6.5  in  dimensionless  form.  For  a  pressure  of  12500  p.s.i. 
the  maximum  deflection  is 

=  0.1380  inches 

which  is  only  50%  of  the  1  inear  deflection.  The  410  and  411  elements  yield 
identical  solutions  for  the  linear  case  as  do  420  and  422.  The  reason  is,  of 
course,  that  the  bending  shape  functions  are  the  same  for  the  different 
element  types  in  the  two  series.  Therefore,  it  was  only  necessary  to  compare 
410  and  420  for  the  linear  case.  Figure  6.6  shows  the  convergence  behavior  in 
terms  of  the  maximum  displacement  for  410  and  420.  It  is  noteworthy  that  the 
presence  of  the  mid-side  rotations  gives  rise  to  a  poorer  result  for  the 
coarsest  mesh.  Table  6.1  shows  the  percent  error  as  a  function  of  the  number 
of  elements. 
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Thus,  for  bending  deformations,  the  420  series  are  clearly  superior  to  the  410 
series  so  far  as  accuracy  and  rate  of  convergence  is  concerned.  However, 
these  results  are  of  limited  significance  in  the  wider  context  of  nonlinear 
shell  analysis.  The  nonlinear  case  brings  in  the  effect  of  membrane 
deformation  as  well  as  bending  and  therefore  gives  a  more  complete  picture  of 
the  overall  convergence  of  the  elements.  Figure  6.7  shows  how  the  nonlinear 
solutions,  obtained  using  the  410  element  (4  elements  and  25  elements), 
compare  with  the  theoretical.  Figure  6.8  shows  the  relative  convergence 
behavior  for  410,  411,  420  and  422.  Table  6.2  shows  the  percent  errors  as  a 
function  of  inesh  size.  The  411  element  appears  to  give  the  best  accuracy  for 
the  coarsest  mesh  in  the  nonlinear  case.  This  appears  somewhat  anomalous  in 
comparison  with  the  rest  of  the  elements  except  for  the  fact  that  411  is  the 
only  one  that  has  two  independent  normal  rotations  at  each  corner  and  thus  has 
an  added  degree  of  membrane  shear  flexibility.  It  is  clear  that  in  this 
problem  it  is  the  refinement  of  the  membrane  shape  function  which  is 
responsible  for  improvements  in  accuracy  for  a  given  mesh  size. 

However,  in  comparing  the  convergence  rates  of  different  elements,  we  must 
beware  of  drawing  premature  conclusions  particularly  for  nonlinear  analysis. 
The  ultimate  goal  of  the  structural  analyst  is  to  produce  the  best  accuracy 
commensurate  with  the  cost  of  performing  the  analysis.  The  cost  for  the 
nonlinear  analysis  is  not  only  a  function  of  the  element  complexity  but  also 
of  the  way  in  which  the  element  affects  the  performance  of  the  solution 
algorithm.  This  is  particularly  relevant  for  STAGSC-1  because  of  the 
automatic  adjustment  of  load  step  size  and  refactoring  operations.  Table  6.3 
provides  details  of  the  17  nonlinear  analyses  performed  for  element 
convergence.  It  is  clear  from  the  table  that  much  depends  on  whether  or  not 
the  step  size  is  halved  durinq  the  analysis.  For  example.  Run  No.  12  using  25 
410  elements  took  less  CPU  time  than  Run  9  using  only  9  elements.  The  reason 
was  that  with  the  more  refined  mesh,  solution  convergence  criteria  were 
satisfied  without  cutting  back  the  load  step.  A  similar  situation  occurred 
for  Runs  19  and  20  (element  420)  except  that  Run  19  did  not  need  to  refactor 
whereas  Run  20  did. 
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Thus,  if  the  total  CPU  time  for  each  run  is  plotted  against  solution  accuracy 
(Fiqure  6.9),  a  somewhat  different  view  of  convergence  (in  the  broader  sense) 
is  obtained.  Element  411  appears  to  be  aqain  the  best  performer,  except  that 
for  higher  accuracy  ( <3%  error)  a  finer  inesh  using  element  410  provides  an 
equally  cost-effective  solution.  Element  420  in  this  context  appears  to  be 
the  least  cost-effective  element. 

6.1.3  ELEMENT  EIGENVALUE  ANALYSIS 

The  specification  of  assumed  displacement  fields  in  terms  of  polynomial  shape 
functions  is  the  almost  universal  method  for  obtaining  the  stiffness 
characteristics  of  finite  elements.  A  subtle  disadvantage  of  this  approach  is 
that  the  fundamental  deformat ional  capability  of  the  element  becomes  almost 
totally  obscured  by  the  algebraic  complexity  of  the  functions.  The  so-called 
"natural-mode"  method  of  formulating  finite  elements  was  developed  by  Argyris 
[28]  and  his  co-workers  prior  to  1965  but  did  not  achieve  widespread  use  as 
the  problem  of  defining  such  modes  for  more  complex  elements  far  out-weighed 
the  advantages.  However,  these  prescribed  natural  modes  of  deformation  are 
desiqned  to  be  orthogonal  and  are  therefore  eigenvectors  of  the  resultant 
stiffness  matrix.  Therefore,  if  an  eiqenmode  analysis  of  any  stiffness  matrix 
(however  it  may  be  generated)  is  performed,  the  so-called  natural  modes  are 
obtained  which  are  often  physically  much  more  revealing  than  the  orioinal 
shape  functions.  Moreover,  as  discussed  by  Gallaqher  [29],  the  complete  set 
of  stiffness  matrix  eigenvectors  must  include  the  set  of  riqid  body  modes 
which  are  identified  by  zero  eigenvalues.  Since  there  can  be,  at  most,  six 
rigid  body  modes  any  number  of  zero  eiqenvalues  greater  than  six  indicates  the 
presence  of  undesired  kinematic  degrees  of  freedom.  This  approach  is 
therefore  also  a  test  for  anomalous  behavior;  this  could  be  inherent  in  the 
form  of  the  shape  functions  which  have  been  assumed,  or  could  be  a  result  of 
the  numerical  procedures  used  for  integrating  the  terms  in  the  stiffness 
matrix . 

Two  possible  methods  of  performinq  the  eigenvalue  analysis  can  be  envisaged. 

A  direct  dynamic  eiqensnlnt inn  could,  in  principle,  be  obtained  for  a  single 
element  by  invoking  the  appropriate  option  in  STAGSC-1.  There  are  two 
objections  to  this  approach,  one  theoretical  and  the  other  practical.  The 
practical  difficulty  is  decisive  c  -  ace  STAGSf-1  is  not  able  to  perform  any 
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analysis  on  a  system  which  is  totally  unsupported.  Thus,  "free-free"  modes  of 
vibration  cannot  be  determined  because  the  solution  procedure  fails  in  matrix 
decomposition.  Therefore,  the  complete  set  of  modes,  which  must  include  rigid 
body  modps,  cannot  be  determined  in  this  way.  The  second  objection  is,  of 
course,  that  the  vibration  mode  analysis  will  obtain  the  eigensolut ion  for  the 
dynamical  matrix  and  not  the  stiffness  matrix  per  se. 

The  second  approach  (the  method  adopted  for  this  work)  is  to  obtain  the  actual 
stiffness  generated  by  the  program  for  a  single  element  and  perform  a  separate 
eigenmode  analysis  on  this  matrix.  STA6SC-1  saves  the  element  stiffness 
matrix  on  TAPE23  for  a  restart  or  on  TAPE8  for  an  eigenvalue  analysis  based  on 
a  nonlinear  stress  state  (either  buckling  or  vibrations).  Since  the  data  on 
TAPF8  is  written  as  a  straightforward  unformatted  write  operation,  the  latter 
method  was  chosen.  The  eigensolution  for  the  stiffness  matrix  was  then 
obtained  by  using  a  specially  written  program  ELMOD  (Appendix).  This  program 
reads  in  the  stiffness  matrix  (written  as  a  lower  triangle)  and  also  a 
solution  vector.  After  the  eigenvalues  and  eigenvectors  have  been  determined, 
the  program  computes  the  strain  energy  associated  with  each  mode  corresponding 
to  the  given  solution  vector.  The  Appendix  provides  a  description  of  ELMOD 
and  a  listing. 

The  usefulness  of  this  method  depends  entirely  on  the  ability  to  associate  the 
individual  terms  of  the  eigenvectors  with  the  corresponding  degrees  of  freedom 
for  the  element.  This  proved  to  be  possible  only  for  the  410  element.  All 
other  elements  were  formed  with  some  degrees  of  freedom  condensed  out  and  the 
identification  of  the  remaining  freedoms  could  not  be  performed. 

Ten  separate  load  cases  were  devised  for  the  single  element  model.  Figure 
6.10  shows  the  model  geometry  and  the  support  conditions  (which  are  sufficient 
only  to  eliminate  rigid  body  freedoms).  The  load  cases  were  chosen  so  as  to 
excite  as  many  different  modes  of  deformation  as  possible  and  also  to  use  as 
many  load  options  as  possible.  Figure  6.11  shows  the  nodal  equilibrium  forces 
printed  out  by  STAGSC-1  for  each  of  the  cases.  Thus,  Cases  1  and  2  are 
identical  in  distribution  (not  in  magnitude)  and  have  displacement  solutions 
which  differ  only  by  a  common  factor.  Cases  3  and  4  are  similar  but  Case  4 
(distributed  edge  load)  produces  self  equilibrating  moments  about  the  normal 
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at  the  two  nodes  on  the  loaded  edge.  Cases  6  and  7  which  are  comparable  cases 
(but  with  the  loading  in  the  Y-direction)  do  not  produce  moments  about  the 
normal.  However,  as  is  clear  from  Figure  6.11,  the  actual  pattern  of  loads 
and  reactions  is  different  because  of  the  support  conditions.  Cases  6  and  7 
yield  identical  displacement  solutions.  Case  5  (uniform  tangential  shear 
loading)  is  very  similar  to  6  and  7  with  respect  to  nodal  force 
distributions.  The  displacements  in  the  X-direction  are  identical  with  6  &  7 
with  only  minor  (but  real)  differences  in  the  Y-direction.  Cases  8  and  9 
(concentrated  corner  moment)  also  produce  identical  displacement  patterns. 

Case  10,  with  a  single  concentrated  moment  about  the  normal  at  the  free  corner 
produced  zero  displacements  and  nodal  equilibrium  forces;  the  obvious 
conclusion  for  Case  10  is  that  STAGSC-1  calculates  no  contributions  to  the 
force  vector  from  moment  components  about  the  normal. 

In  order  to  obtain  the  modal  distribution  of  strain  energy,  the  stiffness 
eigenmodes  were  determined.  Figures  6.12-a  through  6.12-c  provide  qualitative 
sketches  of  the  mode  shapes  calculated.  A  total  of  24  modes  were  obtained 
corresponding  to  the  24  degrees  of  freedom  of  the  stiffness  matrix.  The 
eigenvalues  and  mode  descriptions  are  contained  in  Table  6.4.  The  most 
striking  observation  from  Table  6.4  is  that  there  are  seven  zero  eigenvalue 
modes.  Since  there  can  be  only  six  genuine  rigid  body  modes,  the  existence  of 
an  extra  mode  indicates  the  presence  of  some  spurious  kinematic  freedom.  This 
is  presumably  associated  with  the  inclusion  of  the  rotations  about  the  surface 
normal  at  the  corners  as  separate  degrees  of  freedom.  Inspection  of  the 
eigenvectors  shows  that  for  the  zero  eigenvalue  modes,  the  displacements  and 
rotations  are  all  mutually  consistent  (and  also  with  all  six  rigid  body 
freedoms),  with  the  exception  of  the  rotations  about  the  normal.  These  are  an 
independent  set  and  hence  explain  the  presence  of  the  seventh  zero  mode.  This 
observation  presents  difficulties  in  explaining  the  declared  purpose  of  the 
normal  rotation  [2],  which  is  to  provide  for  cubic  variation  of  in-plane 
displacements  along  an  edge,  i.e.,  if  the  rotations  are  uncoupled  in  the  rigid 
body  modes,  how  are  they  coupled  to  the  deformation  behavior? 

Setting  aside  this  dilemma  for  the  time  being,  attention  will  be  directed  to 
the  deformation  modes.  The  lowest  eigenvalue  (least  stiff)  deformation  mode 
(Mode  8)  is,  somewhat  surprisingly,  purely  membrane.  The  eigenvector  shows 
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that  the  dominant  freedoms  are  the  corner  normal  rotations.  There  is  a 
companion  mode  (21)  in  which  the  translations  are  the  dominant  freedoms  and 
for  which  there  is  therefore  much  greater  overall  stretching.  Mode  9  is  the 
essential  linear  twisting  mode.  Modes  10  through  13  represent  various  types 
of  pure  bending  behavior  with  Mode  10  being  the  anticlastic  bending  mode  and 
Mode  13  spherical  bending.  Modes  14  and  15  are  a  combination  of  bending  and 
twisting.  Modes  16  through  22  are  various  membrane  modes  which  include  pure 
in-plane  shear  (18)  and  pure  membrane  dilation  (22).  The  two  remaining  modes 
are  doubly  antisymmetric  diagonal  bending  and  twisting  modes.  In  Mode  23  the 
diagonals  remain  straight,  while  Mode  24  has  anticlastic  bending  along  the 
diagonals. 

Table  6.5  shows  the  results  of  the  modal  energy  computations.  As  discussed 
previously,  there  are,  in  reality,  only  six  load  cases  which  provide 
distinctly  different  displacement  solution  vectors;  these  are  Cases  1,  3,  4,  6 
and  8. 

Case  1,  as  might  be  expected,  deforms  the  element  almost  entirely  into  a 
linear  twist  mode  (9);  the  only  other  mode  involved  is  24  but  its  energy 
content  is  so  low  that  it  could  be  attributed  to  numerical  round-off.  Load 
Cases  3,  4,  5  and  6,  being  all  in-plane  loadings,  excite  only  membrane  modes. 
The  most  striking  observation  is  that  most  of  the  energy  is  stored  by  Modes  8 
or  16  with  the  other  modes  participating  to  provide  the  asymmetric  response. 
Case  4  is  outstanding  because,  although  the  loading  is  in  the  same  direction 
as  Case  3,  Mode  8  does  not  participate  at  all.  This  may  be  related  to  the 
fact  that  Case  4  is  the  only  one  in  which  non-zero  equilibrium  moments  are 
printed  out  for  the  normal  rotation  directions.  The  response  for  load  Case  8 
is  mainly  provided  by  Modes  9,  10  and  11. 

None  of  the  load  cases  produced  any  significant  response  in  the  diagonal 
bending  plus  twist  modes  (23  &  24)  although  it  is  obvious  that  a  superposition 
of  load  Cases  8  and  9  ought  to  produce  a  response  in  which  Mode  23 
participates  more  strongly.  Accordingly,  the  solution  vectors  for  Cases  8  and 
9  were  superposed  and  the  energy  distribution  calulated.  As  expected.  Mode  23 
now  provides  the  dominant  response. 
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It  is  clear  that  the  only  questionable  result  of  the  eigenvalue  analysis  is 
the  existence  of  the  seventh  zero  eigenvalue  mode  for  element  410.  Eigenvalue 
analyses  were  performed  also  for  elements  411,  420  and  422  despite  the  fact 
that  no  modal  decomposition  of  the  strain  energy  could  be  performed.  None  of 
these  elements  had  more  than  six  zero  eigenvalue  modes.  The  developers 
recommend  that  when  element  410  is  used  at  least  one  normal  rotation  degree  of 
freedom  should  be  constrained. 

6.2  EIGENS0LUTI0N  PERFORMANCE 

Three  particular  aspects  of  performance  were  chosen  for  investigation.  These 
were: 

A.  ability  of  the  algorithm  to  discriminate  between  closely-spaced  or 
multiple  eigenvalues,  and  the  orthonormal ity  of  the  corresponding 
eigenmodes; 

8.  degradation  of  accuracy  with  increasing  mode  number;  and 
C.  convergence  of  the  solution 

Three  models,  of  increasing  complexity,  were  developed  to  examine  these 
questions.  The  first  was  a  simple  cantilever  beam  in  three  dimensions.  For 
equal,  or  slightly  different,  principle  moments  of  inertia  of  the 
cross-section,  multiple  or  closely-spaced  eigenvalues  can  be  easily  obtained. 
The  second  problem  was  a  cantilever  flat  plate  which  is  also  well  documented 
and  therefore  provides  a  further  check  on  accuracy  and  convergence.  Thirdly, 
a  short  cylinder  with  simply  supported  ends  was  chosen  as  an  example  of  a 
shell  structure  with  closely  spaced  modes  at  the  lower  frequencies. 

6.2.1  3-D  CANTILEVER  BEAM 


Figure  6.13  shows  the  geometry  and  finite  element  idealization  of  the  beam. 
Ten  beam  elements  (Element  210)  were  used  along  the  length.  Element  210  has 
cubic  interpolation  for  transverse  displacement  and  linear  interpolation  for 
axial  displacements  and  twist.  The  exact  solution  for  flexural  frequencies 
was  obtained  from  Reference  30  and  is  as  follows: 


09208-868:2 
(33034)  10 


88 


6.1 


f  ■  i£il  [E!/„AlV'? 


where 

nL  is  obtained  as  the  solution  of 


cosh  (nL)  cos  (nL)  +1=0 


6.2 


and 

p  =  density  in  mass  units 
A  =  area  of  cross-section 

The  STAGSC-1  analysis  consisted  of  five  separate  runs.  In  all  cases  only 
lateral  motion  was  permitted.  The  first  two  analyses  confined  vibrations  to 
the  XZ  plane.  The  next  three  permitted  motion  in  both  planes  but  varied  the 
ratio  of  the  principal  moments  of  inertia  in  order  to  observe  the  performance 
of  the  algorithim  in  the  presence  of  closely  spaced  or  multiple  frequencies. 
The  only  difference  between  the  first  two  analyses  was  that  five  modes  were 
determined  in  the  first  case  and  ten  modes  in  the  second.  Table  6.6  contains 
the  frequency  results  for  vibrations  in  the  XZ  plane. 

Comparison  of  the  STAGSC-1  frequencies  with  the  exact  results  shows  that  the 
accuracy  begins  to  degrade  significantly  after  3  or  4  modes.  However,  this  is 
probably  more  a  consequence  of  the  mesh  than  anything  else,  since  there  are 
only  20  active  decrees  of  freedom.  The  results  also  indicated  that  specifying 
S  or  10  modes  for  determination  did  not  affect  the  accuracy  with  which  they 
were  obtained  nor  did  it  affect  the  number  of  iterations  for  convergence  (6  in 
each  case). 

For  the  cases  of  closely  spaced  and  multiple  modes,  results  are  presented  for 
the  runs  made  with  1%  difference  in  cross-section  dimensions  and  equal 
dimensions.  Frequencies  are  tabulated  in  Table  6.7  for  the  first  six  modes  in 
each  case. 

The  determination  of  the  modes  and  frequencies  was  performed  with  accuracy 
equal  to  the  case  where  they  were  not  closely  spaced.  The  only  noticeable 
difference  is  in  the  solutions  for  the  eigenvectors,  where  the  modes  become  a 
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mixture  of  components  in  the  XZ  and  XY  planes.  However,  in  the  case  where 
there  is  only  1%  difference  in  the  frequency  pairs,  the  dominance  of  one  mode 
over  the  companion  is  extremely  strong  (1/107).  For  equal  frequencies,  this 
dominance  is  much  less  but  in  this  situation  no  dominance  is  really  to  be 
expected.  Convergence  of  the  iterative  procedure  was  the  same  in  all  cases. 

The  comparison  with  the  exact  results  indicated  that  STAGSC-1  underestimated 
all  frequencies.  The  exact  solution  is  based  on  classical  beam  theory  with  no 
allowance  for  shear  deformation  or  rotary  inertia.  These  effects  are  not 
included  in  the  beam  elements  so  it  must  be  concluded  that  the  lowering  of  the 
frequency  is  associated  with  the  element  mass  matrix  (which  is  lumped).* 

6.2.2  FLAT  PLATE  CANTILEVER 

This  problem  is  well  documented,  since  theoretical  and  finite  element 
solutions  to  this  problem  are  available  in  References  31  (page  550)  and  26 
(Volume  E,  page  E4,  1-4)  and  it  is  therefore  well  documented.  The  problem 
description  is  contained  in  Figure  6.14.  Table  6.8  shows  a  comparison  between 
the  STAGSC-1  model,  using  element  410,  the  MARC  solution  [26]  using  MARC 
element  4,  the  Zienkiewicz  solution  [31]  using  a  non-conforming  triangular 
element  and  the  classical  solution.  In  each  case  the  mesh  consists  of  two 
square  regions  which  gives  two  elements  for  STAGSC-1  and  MARC  and  four 
triangles  for  the  Zienkiewicz  solutions. 

For  two  elements,  the  STAGSC-1  results  are  very  poor.  This  is  not, 
apparently,  a  failure  of  the  eigensolver  since  much  better  solutions  were 
obtained  using  more  elements. 

Figure  6.15  shows  the  convergence  of  the  first  mode  frequency  for  STAGSC-1  as 
a  function  of  the  number  of  elements.  For  the  purposes  of  comparison  the 
Zinkiewicz  (triangle)  solution  is  plotted  as  a  2-element  mesh  because  the  same 
basic  rectangular  grid  is  used. 

♦Convergence  should  be  from  above.  The  fact  that  the  frequencies  are  under¬ 
estimated  indicates  a  failure  to  satisfy  some  equations  of  mechanics.  It 
may  also  be  the  result  of  numerical  integration  procedures  for  stiffness 
formation. 
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The  eigenvalue  solution  (4  frequencies  requested)  required  4  iterations  for 
the  2  element  model,  6  iterations  for  8  elements  and  7  iterations  for  32 
elements.  Thus,  the  number  of  iterations  required  was  not  strongly  dependent 
on  the  size  of  the  model. 

6.2.3  SIMPLY  SUPPORTED  CYLINDER 

An  analytical  solution  for  the  free  vibrations  of  a  thin-walled  cylinder  was 
first  documented  by  Baron  &  Bleich  [34]  and  their  solutions  were  utilized  by 
other  workers  [32,  33]  for  comparison  with  a  finite  element  solution  and  also 
in  supersonic  flutter  calculations.  A  characteristic  of  the  modal  behavior  of 
a  thin  cylinder  is  that  the  lowest  frequency  is  obtained  for  an  axial 
wavelength  of  2  diameters  {m  =  1)  and  a  circumferential  wavelength  of  */8 
diameters  (i.e.,  8  circumferential  waves;  n  =  8).  For  numbers  of 
circumferential  waves  both  greater  or  less  than  8,  the  frequencies  are 
higher.  Figure  6.16  shows  the  geometry  selected  for  the  STAGSC-1  analysis, 
which  is  the  same  as  that  used  by  Greene,  et.  al.,  [32]  and  Voss  [33].  Since 
the  length  was  chosen  to  be  one  diameter  and  simply  supported  boundary 
conditions  were  imposed  at  each  end,  solutions  can  be  obtained  for  axial 
wavelengths  which  are  equal  to  2D/m,  where  m  is  the  number  of  axial  half  waves 
between  the  ends  and  D  is  the  diameter.  Figure  6.17  shows  how  the  frequencies 
are  distributed  with  respect  to  the  number  of  full  circumferential  waves  n. 

The  lowest  frequency  corresponds  to  a  single  axial  half  wave  (m  =  1),  and 
eight  full  circumferential  waves  (n  =  8).  Attention  was  confined  to  modes 
where  m  =  1  and  therefore  it  was  possible  to  impose  symmetry  boundary 
conditions  at  the  mid  axial  section.  Some  initial  runs  were  made  in  order  to 
establish  a  reasonable  mesh  for  the  analysis.  Greene,  et.  al.,  [32]  analysed 
the  problem  using  a  mesh  based  on  a  quarter  wave  model,  i.e.,  the  model 
spanned  11  1/4°  circumferentially  and  half  the  length,  which  is  1/4  wave  in 
both  directions  for  the  fundamental  mode  (m  =  1,  n  =  8,  Figure  6.13).  For  the 
present  study,  two  1/4  wave  meshes  were  set  up,  one  with  a  4  x  4  mesh  (16 
elements)  and  one  with  a  5  x  5  mesh.  The  fundamental  frequencies  obtained  are 
given  in  Table  6.9.  These  results  indicate  a  wel 1 -converged  solution,  at 
least  for  the  first  few  modes.  However,  the  11  1/4°  model  is  only  suitable 
for  the  fundamental  (or  appropriate  multiples)  because  of  the  boundary 
conditions  imposed  on  the  axial  edges.  In  order  to  establish  a  mesh  suitable 
for  higher  modes,  a  1/8  cylinder  model  was  set  up. 
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Referring  to  Figure  6.16,  the  1/8  model  spanned  the  region  ABCD;  symmetric 
boundary  conditions  were  imposed  on  edges  AB,  BC  and  AD.  Table  6.10  presents 
results  for  four  different  meshes. 

The  results  show  that  in  order  to  capture  the  lowest  modes  accurately,  at 
least  16  circumferential  elements  were  required  (i.e.,  2  per  quarter  wave)  and 
4  axial  (4  per  quarter  wave). 

The  larger  number  per  quarter  wave  in  the  axial  direction  is  probably  due  to 
the  poor  aspect  ratio  obtained  with  the  16  x  2  mesh.  This  seems  to  be 
substantiated  by  comparing  the  16  x  4  and  the  16  x  8  meshes.  Figure  6.18 
shows  the  corresponding  mode  shapes  and  emphasizes  the  fact  that  the  best 
accuracy  is  obtained  for  Mode  5  (n  =  6)  where  the  best  definition  of  the  mode 
shape  is  obtained.  Thus,  in  order  to  develop  reasonable  accuracy  for  a  full 
circumference  model,  it  became  clear  that  64  circumferential  elements  would  be 
necessary  and  at  least  2  axial  elements  for  a  half  cylinder  model  (using 
symmetry  about  the  mid-length). 

Using  the  1/2  cylinder  model,  some  apparently  anomalous  results  were 
obtained.  The  first  two  modes  for  the  cylinder  correspond  to  n  =  8  and  n  =  7 
circumferential  waves  respectively.  STAGSC-1  was  executed  with  the  following 
parameters  specified 

NEIG  =  4,  SHIFT  =  E IGA  =  EIGB  =  0 

which  instructs  the  program  to  determine  the  four  lowest  eigenvalues  using  a 
zero  frequency  shift.  The  frequencies  obtained  are  shown  in  Table  6.11. 

The  modes  that  were  isolated  by  the  program  consisted  of  two  pairs  each  with  8 
and  7  circumferential  waves  respectively,  but  with  slightly  different 
frequencies.  Such  pairs  do  not  exist  according  to  the  classical  analysis. 
Figure  6.19  through  6.22  show  the  mode  shapes  in  terms  of  the  normal 
displacements  plotted  around  the  circumference.  It  seems  clear  that  modes  1 
and  3  (as  numbered  by  STAGSC-1)  correspond  to  the  true  modes  with  8  and  7 
circumferential  waves  whereas  modes  2  and  4  are  anomalous  due  to  their 
irregularity  in  circumferential  distribution.  The  eigensolver  obtained 
converged  eigenvalues  after  a  total  of  16  iterations,  having  automatically 
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selected  a  frequency  shift  of  78.53  Hz  after  8  iterations.  Thus,  the 
alqor i thin  was  quite  capable  of  discrimination  between  the  closely-spaced 
frequencies  in  each  pair.  It  seems  likely,  therefore,  that  the  anomalous 
modes  are  a  result  of  the  finite  element  discretization.  Since  the  element 
type  and  mesh  size  are  the  same  as  for  the  1/8  cylinder  model  (which  did  not 
give  rise  to  these  paired  modes)  the  most  likely  source  of  the  anomaly  is  the 
juncture  of  the  cylindrical  shell  unit.  This,  in  effect,  introduces  an  axial 
"seam"  in  the  structure  which  probably  results  in  some  slight  asymmetry. 

If  this  is  indeed  the  situation,  then  it  could  be  argued  that  the  eigensolver 
is  highly  sensitive  and  discriminatory  in  identifying  modes  which  are  so 
closely  similar.  On  the  other  hand,  the  more  practical  conclusion  is  that 
such  powers  of  discrimination  may  be  a  considerable  nuisance  where,  in  an 
unknown  situation,  the  separation  of  real  and  spurious  modes  may  not  be  so 
easy. 

In  order  to  obtain  further  modes,  use  was  made  of  the  feature  whereby  the 
frequencies  in  the  vicinity  of  a  given  range  can  be  determined  (EIGA,  EIGB). 
The  complete  set  of  results  is  presented  in  Table  6.12. 

6.2.4  CONCLUSIONS 

The  evaluation  exercises  performed  permit  some  overall  conclusions  to  be  drawn 
with  respect  to  the  performance  of  the  eigensolver.  Three  aspects  were 
investigated  and  these  can  be  summarized  as  follows: 

A.  Discrimination  between  closely  spaced  or  multiple  modes  is  very  good 
for  problems  of  widely  varying  complexity. 

B.  Degradation  of  accuracy  with  increasing  mode  number  was  observed,  but 
the  underlying  cause  was  probably  associated  more  with  the  adequacy 
of  the  mesh  and  the  elements  themselves,  rather  than  the  eigensolver. 

C.  Solution  convergence  behavior  for  the  problems  investigated  was 
satisfactory. 

93 
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Overall  conclusions  based  on  this  study  are  therefore  that  the  eigensolver  in 
STAGSC-1  is  highly  discriminatory,  accurate  and  efficient.  In  fact,  the 
quality  of  the  solutions  obtained  was  probably  mostly  dependent  on  the 
elements  themselves  rather  than  the  eigensolver. 

6.3  TRANSIFNT  INTEGRATION  PERFORMANCE 

STAGSC-1  is  particularly  well  equipped  with  transient  integration  capability 
having  one  explicit  and  four  implicit  operators.  The  operators  themselves  are 
described  in  Section  4.8,  and  it  is  the  purpose  of  this  section  to  investigate 
their  relative  performance.  Much  has  been  written  about  the  accuracy, 
stability,  damping  and  so  on  of  the  multitude  of  numerical  operators  which 
have  been  developed  for  the  integration  of  ordinary  differential  equations. 
Most  of  the  research,  however,  has  dealt  with  their  application  to  linear 
problems,  presumably  because  the  mathematical  proofs  involved  can  be  more 
readily  derived.  From  the  practical  point  of  view,  use  of  these  algorithms  in 
a  program  such  as  STAGSC-1,  more  often  than  not,  will  be  for  the  integration 
of  a  non-linear  set  of  equations.  The  main  thrust  of  the  present  study  was 
therefore  directed  towards  studying  the  performance  of  the  operators  when 
applied  to  a  non-linear  problem.  The  properties  which  are  relevant  are  the 
same  as  for  a  linear  problem,  i.e.,  accuracy,  stability,  artificial  damping 
and  frequency  distortion  (dispersion). 

Two  problems  were  selected  for  the  investigation.  The  first  was  a  benchmark 
to  establish  the  general  validity  of  the  transient  integration  in  STAGSC-1, 
and  was  chosen  to  be  the  linear  elastic  response  of  a  thin-canti lever  flat 
plate  subjected  to  a  triangular  pressure  pulse.  This  example  is  a 
demonstration  problem  used  for  the  MARC  program  [26]  and  a  comparison  solution 
was  therefore  available.  The  second  problem  was  the  non-linear  response  of  a 
circular  ring  segment  subjected  to  an  impulsive  pressure  loading.  The 
solution  to  this  problem  has  been  reported  by  Stricklin,  et.  al.,  [37]  with 
comparisons  between  analysis  and  experiment.  All  of  the  candidate  integration 
operators  were  tested  using  this  problem. 
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6.3.1  FLAT  PLATE  MODEL 


The  basic  model  and  dimensions  were  the  same  as  were  used  in  the  work  on  the 
eigensolver  performance  and  details  may  be  found  in  Figure  6.14.  The  pressure 
time  history  is  shown  in  Figure  6.23  and  is  a  triangular  shaped  pulse  with  a 
peak  value  of  100  lb/in  applied  uniformly  over  the  whole  plate.  The  total 
duration  of  the  pulse  was  0.04  milliseconds. 

This  problem  was  analyzed  using  the  explicit  and  the  implicit  trapezoidal 
methods  (Newmark  -  a)  in  STAGSC-1.  The  model  used  element  410  in  a  mesh 
with  four  elements  lengthwise  and  two  across  the  width.  The  bench  mark 
comparison  was  obtained  using  the  MARC  program.  The  MARC  model  was  based  on 
element  4  and  a  mesh  consisting  of  two  elements  lengthwise. 

The  explicit  method  requires  the  selection  of  a  time  step  which  is  below  the 
stability  limit.  This  may  be  estimated  in  several  different  ways  which  may  be 
summarized  as  follows: 


where  <«max  is  the  highest  circular  frequency  which  is  inherent  in  the 
finite  element  model. 
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where 


C?  =  E/p  ( 1  -  v2 ) ; 

=  G/p 

au  (<_  Aa)  are  grid  spacings 
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h  -  interval  in  finite  difference  formulation;  iri  this  case  the 
assumption  is  made  that  Aa  =  2h 

(see  [3],  page  6-33) 


At 


Act 

<  — 
-  C 


6.4 


where 


?  ? 

C  =  E/p  ( 1  —  v  )  =  speed  of  sound 

llsinq  the  material  properties  for  the  flat  plate  model  and  the  spacing 
Aa  =  ab  =  0.5  inch  (4x2  mesh)  the  following  values  of  At  were  obtained 
and  are  tabulated  below; 

The  table  shows  that  a  At  of  less  than  2  x  10"6  seconds  is  required  for  a 
stable  solution.  Accordingly  a  time  step  of  1  x  10"6  seconds  was  chosen  but 
the  solution  diverged  after  6  time  steps.  A  tenfold  reduction  in  time  step 
size  to  0.1  x  10  6  seconds  permitted  the  solution  to  progress  further  but 
divergence  still  occurred  after  13  time  steps.  The  most  probable  source  of 
this  discrepancy  in  the  assumption  used  in  the  critical  time  step  formula 
6.3.  For  high-order  finite  elements,  the  values  of  h  (finite  difference 
interval)  or  ao  are  difficult  to  estimate,  and  these  estimates  can  be 
seriously  affected  by  inconsistent  mass-lumping  combined  with  high-order 
deformation  patterns.  As  pointed  out  by  Krieg  and  Key  [35],  explicit 
operators  do  not  work  effectively  with  high-order  elements  for  this  reason. 
However,  the  fact  remained  that  the  time  step  size  clearly  had  to  be  much 
smaller,  probably  by  at  least  an  order  of  magnitude,  in  order  that  a  stable 
solution  be  obtained.  The  effort  was  therefore  discontinued  at  this  point  and 
attention  directed  towards  the  implicit  trapezoidal  operator. 

The  MARC  solution  referred  to  earlier  had  used  a  time  step  of  20  x  10-^ 
seconds.  Accordingly,  for  the  STAGSC-1  solution  values  of  10  x  10”^ 
seconds,  20  x  10  ^  seconds  and  100  x  10  ^  seconds  were  used  and  the 
solutions  compared  with  the  MARC  response.  Solutions  were  obtained  out  to  a 
total  time  of  1000  x  10  ^  seconds.  Figures  6.24  and  6.25  show  the  tip 
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displacement  and  tip  velocity  as  a  function  of  time  for  the  STAGSC-1  and  MARC 
solutions.  First,  there  was  no  significant  difference  in  the  STAGSC-1  tip 
displacement  time  histories  for  at  =  10  x  10”^  seconds  and  at  =  20  x 
10'6  seconds  so  results  are  presented  only  for  at  =  10  x  10  ^  seconds. 

Figure  6.24  shows  that  there  is  excellent  agreement  between  the  MARC  and 
STAGSC-1  solutions.  At  first  sight  there  appear  to  be  some  dispersive  affects 
in  the  STAGSC-1  solution  as  compared  with  MARC  since  the  fundamental  response 
frequency  is  about  14%  less  than  the  MARC  result.  However,  it  is  more  likely 
that  the  "dispersion"  is  merely  underprediction  of  the  frequencies  by  the 
STAGSC-1  model  (Element  410)  as  previously  demonstrated.  The  peak-to-peak 
amplitude  of  the  STAGSC-1  solution  is  about  the  same  as  that  predicted  by 
MARC.  The  solution  obtained  using  a  time  step  of  100  x  10"^  seconds  is 
clearly  rather  coarse  and  does  not  give  any  resolution  of  the  higher  frequency 
components  as  might  be  expected.  However,  the  solution  appeared  to  be  on  the 
point  of  diverging  at  the  final  time  step.  Figure  6.25  shows  how  the 
responses  in  terms  of  tip  velocity  compare.  The  previously  observed  frequency 
distortion  is  rather  more  pronounced,  particularly  with  regard  to  the  higher 
frequency  components.  Thus,  it  may  be  inferred  from  this  benchmark  case  that 
correct  and  reasonably  accurate  results  car*  be  obtained  using  the  STAGSC-1 
trapezoidal  operator  for  a  linear  response  analysis.  The  failure  to  do  so 
with  the  explicit  method  could  probably  be  resolved  by  using  a  smaller  time 
step,  but  this  was  not  demonstrated. 

6.3.2  IMPULSIVELY  LOADED  RING 

Stricklin,  et.  al.,  [37]  performed  a  large  deflection  elastic-plastic 
transient  dynamic  analysis  of  an  impulsively  loaded  ring  using  the  program 
DYNAPLAS.  The  ring  problem  was  initially  posed  and  analyzed  by  Wu  &  Witmer 
[38]  who  compared  their  results  with  experimental  data.  The  geometry  of  the 
ring  is  shown  in  Figure  6.26.  The  initial  conditions  for  the  problem  are  zero 
displacements  and  a  radially  inward  initial  velocity  of  4862  inch/second 
imposed  over  a  central  sector  of  120°.  The  STAGSC-1  finite  element  model 
assumed  symmetry  about  the  plane  through  the  middle  of  the  ring  (a  =  0°)  and 
utilized  21  equal  elements  between  e  =  0°  and  e  =  157.5°.  The  element 
used  was  element  410. 
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The  original  intent  was  to  perform  full  elastic-plastic  large  deflection 
transient  analyses  for  each  operator  and  for  various  time  steps.  However, 
difficulties  were  encountered  with  the  plasticity  solution  and  this 
comprehensive  approach  was  discarded  in  favor  of  an  elastic  large  displacement 
solution.  This  was  considered  to  be  sufficiently  nonlinear  to  test  the 
integration  operators  effectively.  Table  6.14  gives  a  summary  of  the 
individual  analyses  performed. 

Runs  10  through  11  were  all  full  non-linear  elastic-plastic  analyses  using  the 
Park  integration  operator. 

A.  Park  Operator 

A  time  step  of  2  x  10  ^  seconds  was  chosen  for  the  first  run  on  the 
assumption  that  this  would  be  sufficiently  small  to  permit  accurate 
handling  of  the  plasticity.  The  results  indicated  that  the  first  40 
steps  required  plastic  sub-iterations  at  each  step  but  thereafter  the 
need  became  only  occasional.  Solution  convergence  was  obtained  with 
only  a  single  iteration  per  time  step  after  the  first  100  steps.  The 
solution  reached  a  total  time  of  620  x  10"6  seconds  when  the 
execution  time  limit  was  reached.  Accordingly,  the  analysis  was 
repeated  using  the  automatic  time  step  control  with  an  initial  step 
of  2  x  10  6  seconds  imposed  for  the  first  100  steps.  Subsequently, 
the  step  was  doubled  twice  during  the  next  16  steps.  A  third 
doubling  was  attempted  but  the  solution  failed  because  the 
determinant  changed  sign.  It  appears  that  in  the  automatic  mode 
STAGSC-1  at  present  is  able  only  to  increase  the  time  step  which  is 
obviously  unsatisfactory  for  non-linear  analysis  in  general  where  the 
changing  character  of  the  nonlinearities  as  the  solution  progresses 
may  require  successive  increases  or  decreases  in  the  step  size. 

For  the  third  run  (Run  12),  a  fixed  step  of  2  x  10‘6  seconds  was 
used  and  a  restart  file  (TAPE22)  was  saved.  The  solution  was 
obtained  up  to  a  time  of  1002  x  10'^  seconds. 
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The  fourth  run  was  an  attempted  restart  which  failed  giving  the 

message;  "Restart  tape  does  not  contain  plasticity  data  .  . 

notwithstanding  the  fact  that  the  previous  run  generating  the  file 

confirmed  that  the  plasticity  data  had  been  written  to  tape.  At  this 

point,  the  solution  was  examined  in  detail  and  compared  with  the 

published  results.  Figure  6.27  shows  the  time  history  of  the 

deflection  at  the  center  of  the  ring  up  to  0.001  seconds.  The 

solution  compares  reasonably  well  with  the  DYNAPLAS  solution  [37]  but 

shows  a  qrowing  difference.  The  circumferential  strain  history 

(Fiqure  6.28)  at  the  center  of  the  ring  (outer  surface)  also  shows 

good  agreement.  Perhaps  the  most  encouraging  comparison  is  with  the 

experimental  data  presented  by  Wu  and  Witmer  [38].  The  deformed 

shape  at  0.0008  seconds  as  calculated  by  STAGSC-1  is  compared  with 

the  experimental  results  at  0.00079  seconds  in  Figure  6.29. 

Nevertheless,  it  was  concluded  that  for  times  greater  than  this  the 

solution  was  flawed.  Figure  6.30  compares  the  distribution  of 

bending  rotation  at  0.001  seconds  for  STAGSC-1  and  DYNAPLAS. 

STAGSC-1  shows  large  discontinuities  in  rotation  at  three  locations 

(b  =  52.5°,  75°  and  150°)  which  clearly  are  unacceptable.  This 

appeared  to  be  traceable  to  the  plasticity  solution  since  the 

stresses  did  not  lie  on  the  stress-strain  curve.  Since  the  emphasis 

was  intended  to  be  on  the  integration  operators,  it  was  decided  at 

this  point  to  avoid  the  plasticity  problems  and  to  perform  the 

evaluation  using  a  purely  elastic  large  displacement  model.  Thus, 

runs  14  through  25  (Table  6.14)  were  confined  to  the  elastic  regime. 

Using  the  same  time  step  as  for  the  plasticity  solution,  response  was 

obtained  up  to  the  maximum  time  specified  of  0.002  seconds.  The 

solution  was  considered  to  be  quite  accurate  since  only  one  or  two 

iterations  per  time  step  were  required.  For  the  next  run,  the  time 

-6 

step  was  increased  by  a  factor  of  ten  (20  x  10*  seconds).  This 
produced  an  interesting  result;  after  the  first  step,  the  time 
increment  was  automatical ly  halved  and  the  stiffness  matrix 
refactored,  in  contrast  with  run  11  where  automatic  time  step  control 
was  used  and  the  time  step  was  not  reduced  when  difficulties  were 
encountered.  Since  a  small  time  limit  had  been  imposed,  the  run  was 
repeated  with  a  time  step  of  10  x  10"6  seconds.  This  solution 
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refactored  periodically  and  required  from  2  to  8  iterations  per  time 
step.  Figure  6.31  shows  the  responses  obtained  for  the  two  time 
steps.  Increasing  the  time  step  eliminated  some  of  the  high 
frequencies  and  increased  the  dispersion  (about  5%). 

B.  Trapezoidal  Operator 

Time  steps  of  2  x  10"^  seconds  and  10  x  10  ^  seconds  were  again 
chosen.  Figure  6.32  shows  the  responses  obtained.  There  appears  to 
be  very  little  to  distinguish  the  two  solutions;  dispersion  and 
damping  are  both  unaffected  by  increasing  the  time  step.  Comparison 
with  Figure  6.31  shows  that  the  results  are  almost  identical  with 
those  obtained  using  the  Park  operator. 

C.  Gear  2nd  and  3rd  Order  Operators 

The  Gear  2nd  order  operator  shows  about  an  8.4%  dispersive  increase 
in  the  fundamental  period  when  the  time  step  is  increased  (Figure 
6.30).  In  addition,  the  damping  of  the  higher  frequencies  is  very 
marked  in  comparison  with  the  other  operators.  The  Gear  3rd  order 
operator  became  unstable  at  about  .00035  sec.  using  the  smaller  time 
step.  No  further  runs  were  attempted. 

0.  Explicit  Operator 

Time  steps  of  1  x  10"^  seconds,  5  x  10  ^  seconds  and  1  x  10 
seconds  were  attempted.  The  two  larger  steps  both  led  to  divergence 
after  a  small  number  of  steps  (Figure  33).  The  run  using  a  step  of  1 
x  10"7  seconds  did  not  diverge  but  completed  only  1710  steps  (t  = 
.00017  seconds)  when  the  time  limit  for  the  run  was  reached.  Since 
this  was  clearly  extremely  inefficient  compared  with  the  other 
operators,  no  further  runs  were  attempted.  The  details  of  the 
computer  run  times  and  resource  usage  are  presented  in  Table  6.15. 

Examination  of  the  execution  times  presented  in  Table  6.12  shows  that  for  the 
implicit  operators  there  are  no  obvious  grounds  for  selecting  one  over  another 
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on  the  basis  of  computational  efficiency.  Comparison  of  the  Park  and 
trapezoidal  methods  (Runs  14  and  18)  shows  that  when  the  numbers  of  iterations 
are  the  same  and  there  is  no  refactoring  then  the  execution  times  and  computer 
resource  usage  are  also  the  same.  Where  differences  do  occur  (e.g., 
trapezoidal  and  Gear's  2nd  Order  -  Runs  17  and  19)  they  can  be  ascribed  to 
more  refactoring  in  the  case  of  the  one  method  (trapezoidal)  compared  with  the 
other.  Other  things  being  equal,  it  may  be  said  that  this  also  reflects 
properties  of  the  integration  operator,  but  such  a  conclusion  would  need  to  be 
more  firmly  based  than  is  possible  with  the  present  evidence.  Such  effects  as 
problem  dependence  would  need  to  be  investigated. 

Perhaps  the  most  surprising  result  is  the  very  poor  performance  of  the 

explicit  method  which  should  be  far  superior  on  the  basis  of  resources  used 

per  time  step.  This  is  clearly  not  the  case  given  the  present  problem  and  the 

low  value  of  the  initial  time  step  (0.1  x  10"6  seconds  <  at  <  0.5  x 
“6 

10  seconds)  makes  it  hopelessly  uneconomical.  For  a  large  problem 
however,  where  the  cost  of  refactoring  will  be  much  greater,  the  explicit 
method  should  become  competitive. 

6.3.3  CONCLUSIONS 

Overall,  the  STAGSC-1  integration  operators  performed  satisfactorily  for  both 
linear  and  nonlinear  transient  response  problems.  For  the  relatively  small 
test  problems  the  explicit  method  did  not  show  to  advantage  because  the  size 
of  the  critical  time  step  was  so  small  and  the  number  of  steps  became 
excessive. 

Of  the  implicit  methods,  the  trapezoidal  operator  was  the  most  effective  and 
showed  the  smallest  dispersion  when  the  time  step  was  increased.  The  Gear  2nd 
and  3rd  order  operators  were  the  least  effective,  the  3rd  order  method 
becominq  unstable  even  for  the  smaller  time  step. 

6.4  NONLINFAR  COLLAPSF  ANALYSIS 

This  section  presents  results  obtained  from  a  number  of  elastic  nonlinear 
collapse  analyses  of  point  loaded  cylindrical  shells.  The  analyses  were 
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performed  more  or  less  on  a  production  basis  as  part  of  an  investigation  of 
the  buckling  of  a  nuclear  reactor  containment  vessel  due  to  localized  loads, 
and  not  as  part  of  the  STAGSC-1  evaluation  per  se.  However,  because  the 
analyses  are  for  the  less  frequently  considered  case  of  displacement-control , 
and  because  of  experience  gained  with  certain  input  strategy  parameters,  it  is 
fruitful  to  discuss  these  aspects  of  the  analyses  as  they  pertain  to  the 
STAGSC-1  evaluation.  In  this  respect,  it  should  be  pointed  out  that  the 
analyses  might  very  well  have  been  performed  differently,  and  certainly  for 
simpler  and  more  economical  problems,  if  the  objective  at  the  outset  had  been 
an  evaluation  of  the  nonlinear  collapse  capability  of  STAGSC-1. 

The  problem  of  interest  for  the  containment  vessel  buckling  investigation  is 
collapse  of  the  so-called  "poked  cylinder",  which  is  a  cylindrical  shell 
subjected  to  an  inward-directed  normal  point  force  applied  at  midlength.  Two 
related  and  more  standard  problems  were  first  run  on  STAGSC-1  to  gain 
familiarity  with  the  program,  and  to  provide  check  cases.  These  problems  are 
the  point  loaded  Venetian  blind  and  the  pinched  cylinder.  Dimensions  and 
material  properties  for  the  three  elastic  shell  problems  are  given  in  Table 
6.16  along  with  the  STAGSC-1  predictions  of  the  collapse  loads.  Table  6.16 
also  serves  to  define  some  of  the  notation. 

6.4.1  ANALYSIS  CONSIDERATIONS 

Either  displacement  or  load  can  be  specified  as  the  controlled  variable  in 
STAGSC-1  input.  For  reasons  of  economy,  and  to  obtain  the  post-collapse 
behavior,  displacement-controlled  analyses  were  performed  for  each  of  the 
three  shell  problems.  Specifically,  the  radial  displacement  wq  under  the 
point  load  P  is  specified  in  steps,  and  the  corresponding  value  of  P  required 
to  impose  the  specified  displacement  is  determined  from  a  printout  of 
"equilibrium  forces."  The  nonlinear  collapse  load  is  the  value  of  the  load  at 
the  maximum  point  (limit  point,  dP/dwQ=0)  of  the  load-displacement(P-wo) 
curve  obtained  from  the  analysis.  This  load,  which  is  also  called  the 
snap-through  load  or  nonlinear  buckling  load  in  the  literature,  should  not  be 
confused  with  the  bifurcation  buckling  load  obtained  from  an  eigenvalue  type 
of  analysis.  For  comparison,  a  load-controlled  collapse  analysis  was  also 
performed  for  one  of  the  problems.  In  this  case,  the  collapse  load  is  taken 
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to  be  the  value  of  P  for  which  the  finite  element  computations  failed  to 
converge  within  the  specified  convergence  criterion  (DELX),  specified 
permissible  number  of  cuts  in  the  load  step  (NCUT),  and  the  number  of  times 
the  stiffness  matrix  can  be  refactored  (NEWT).  Because  of  symmetry 
conditions,  each  shell  problem  is  analyzed  as  a  cylindrical  panel  whose  axial 
length  is  half  the  actual  length  of  the  shell,  and  whose  circumferential  angle 
is  half  the  actual  opening  angle  of  the  Venetian  blind,  and  90°  and  180°  for 
the  pinched  and  poked  cylinders,  respectively.  Classical  simple  support 
(diaphragm)  boundary  conditions  are  imposed  on  the  curved  edges  of  all  three 
panels,  and  free  edge  conditions  are  specified  on  the  remaining  straight  side 
of  the  Venetian  blind  panel.  The  quadrilateral  plate  element  411  with  2x2 
Gaussian  integration  points*  is  used  in  the  analyses.  Table  6.17  gives  the 
number  of  mesh  points  along  with  computer  costs  and  other  details  of  the 
analyses. 

Analyses  of  the  Venetian  blind  and  pinched  cylinder  were  performed  in  a 
straightforward  manner,  with  no  decisions  required  for  selection  of  the 
strategy  parameters  beyond  initial  specification  of  NCUT  and  NEWT. 
Specifically,  it  was  not  necessary  to  tighten  the  convergence  criterion  DELX 

_3 

from  its  default  value  of  10  ,  or  to  force  refactoring  at  each  of  the 

higher  load  steps.  In  contrast,  the  poked  cylinder  analysis  required  that 
greater  attention  be  paid  to  these  strategy  parameters,  with  DELX  being 

_d  _c 

reduced  from  10  to  10  and  10  ,  and  with  refactoring  forced  at 

consecutive  higher  load  steps  (see  Table  6.17),  as  is  discussed  more 
completely  later. 

6.4.2  VENETIAN  BLIND  &  PINCHED  CYLINDER  RESULTS 

A.  Comparison  With  Previous  Solutions 

Table  6.18  shows  that  linear  STAGSC-1  results  appropriate  to  the 
meshes  of  Table  6.17  agree  to  within  10%  with  previous  solutions 
[41,42]  for  wQEt/P,  a  dimensionless  flexibility  coefficient.  Use 


+Use  of  2x2  integration  points  was  recommended  by  the  developers  of  STAGSC-1. 
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of  a  finer  mesh  is  seen  to  reduce  the  discrepancy  in  the  pinched 
cylinder  results  from  10  to  2%.  For  reasons  of  economy,  the  finer 
mesh  was  not  used  in  the  nonlinear  collapse  analysis  of  the  pinched 
cylinder. 

With  regard  to  the  comparison  of  nonlinear  results  for  the  Venetian 
blind  and  pinched  cylinder,  it  is  noted  that  the  present  results  are 
based  on  the  same  meshes  used  in  the  STAGSA  finite  difference 
analyses  performed  by  Brogan  and  Almroth  [41].  The  above  comparison 
of  linear  results  indicates  that  these  meshes  are  not  sufficiently 
fine  to  render  completely  converged  solutions  for  wQ.  Furthermore, 
the  rate  of  convergence  with  mesh  refinement  is  likely  to  be 
different  for  the  STAGSA  and  STAGSC-1  solutions.  Consequently,  the 
STAGSC-1  results  for  the  Venetian  blind  (Figure  6.34)  and  pinched 
cylinder  (Figure  6.35)  agree  as  well  as  could  be  expected  with  the 
STAGSA  solutions,  with  respective  collapse  loads  Pc  for  the 
Venetian  blind  differing  by  about  5%  (Table  6.16).  A  similar 
comparison  of  collapse  loads  for  the  pinched  cylinder  could  not  be 
made  because  the  STAGSA  solution  in  [41]  was  terminated  before 
collapse. 

B.  Nonlinear  Behavior  of  Venetian  Blind 

For  comparison,  Figure  6.34  shows  load-displacement  plots  obtained 
from  both  the  displacement-controlled  (DC)  and  load-controlled  (LC) 
analyses.  The  following  observations  and  comments  stem  from  an 
examination  of  this  figure  and  Tables  6.16  and  6.17: 

(1)  The  collapse  loads  predicted  by  the  DC  and  LC  analyses  agree  to 
three  significant  figures. 

(2)  The  LC  analysis  required  over  8  times  as  many  steps  to  reach  the 
limit  point  and  is  6  times  as  expensive  up  to  this  point. 

(3)  The  relatively  large  number  of  load  steps  provides  a  very  smooth 
LC  curve  at  all  load  levels.  The  DC  curve  drifts  away  from  the 
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more  accurate  LC  curve  in  a  local  region  and  then  properly 
returns  to  it  upon  a  ref actor ization  at  a  later  displacement 
step  (w  =  .02").  Since  the  same  (default)  value  of  the 
convergence  criterion  is  used  in  both  analyses  (see  Table  6.17), 
and  since  there  are  fewer  steps  in  the  DC  analysis,  it  appears 
that  satisfaction  of  convergence  is  so  much  easier  in  a  DC 
analysis  that  it  sometimes  accepts  too  large  a  displacement 
step.  This  apparently  causes  a  "drift"  in  the  DC  solution  from 
the  "true"  solution  (appropriate  to  the  chosen  mesh)  that  is 
eventually  corrected  by  the  Newton-Raphson  solution  procedure. 
Such  a  local  drift  Is  generally  of  little  consequence. 

(4)  STAGSC-1  has  the  advantage  and  convenience  of  automatically 
reducing  the  load  or  displacement  step  (as  permitted  by  NCUT) 
when  convergence  difficulties  are  encountered.  Thus,  it  is 
easier  for  STAGSC-1  to  home-in  on  the  collapse  load  in  a  LC 
analysis  than  it  is  for  other  programs,  such  as  MARC  [26],  for 
example,  in  which  input  values  of  the  load  steps  cannot  be 
reduced  internally  by  the  program  to  reflect  the  increasing 
ill-conditioning  of  the  stiffness  matrix  as  P*PC-  Such 
programs  generally  require  more  restarts  to  obtain  the  collapse 
load  to  the  same  degree  of  accuracy. 

(5)  The  post-collapse  curve  obtained  from  the  DC  analysis  shows  that 
the  point  loaded  Venetian  blind  is  imperfection-insensitive  in 
the  sense  that  collapse  loads  for  perfect  and  slightly  imperfect 
Venetian  blinds  are  expected  to  differ  by  a  small  amount.  Of 
course,  the  post-collapse  curve  cannot  be  determined  from  a  LC 
analysis  unless  special  procedures  are  used. 

C.  Nonlinear  Behavior  of  Pinched  Cylinder 

The  sketch  of  the  12x8  variable  mesh  given  in  [41]  was  scaled  to 
provide  the  axial  and  hoop  intervals  given  in  Table  6.19.  The 
corresponding  STAGSC-1  solution  in  shown  by  the  upper 
load-displacement  plot  of  Figure  6.35.  Inspection  of  the  figure 
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reveals  that  collapse  of  the  pinched  cylinder  is  not  catastrophic . 

The  load  drops  off  by  only  1%  from  its  value  at  collapse  to  its 
minimum  post-collapse  value,  while  the  corresponding  displacement 
increases  by  just  4%  before  stiffening  behavior  takes  place.  This 
portion  of  the  P-wq  plot  in  the  vicinity  of  the  maximum  and  minimum 
points  was  routinely  traced  out  by  STAGSC-1,  as  was  the  entire  P-w 

o 

curve. 

These  results  are  useful  in  a  general  sense  because  this  is  the  first  time 
that  collapse  has  been  predicted  for  a  pinched  cylinder.  Hence  collapse  of 
the  considerably  thinner  poked  cylinder  (R/t  =  638  versus  100}  was  seen  to  be 
a  definite  possibility,  and  was  therefore  investigated.  However,  specific 
numerical  values  and  details  of  the  pinched  cylinder  results  may  be  inaccurate 
because  (1)  the  mesh  is  not  fine  enough  to  provide  converged  values  of  the 
radial  displacement  and  smooth  hoop  spatial  variations  for  stresses,  (2) 
computed  rotations  (^25°  at  collapse)  exceed  the  range  of  permissible  values 
for  valid  application  of  the  theory  employed  in  STAGSC-1,  and  (3)  computed 
stresses  at  collapse  for  this  relatively  thick  cylinder  indicate  that 
plasticity  effects  (not  considered  here)  are  significant  even  at  points  away 
from  the  load. 

6.4.3  ANALYSIS  OF  THE  POKED  CYLINDER 

A.  Mesh  Considerations 

It  was  decided  to  use  the  finest  mesh  that  would  not  be  unduly 
expensive  and  that  would  not  exceed  the  core  storage  limitation  that 
the  total  number  of  nodes  mxn  could  not  exceed  528  (on  the 
Westinghouse  CDC7600),  which  was  determined  by  trial  and  error  from 
linear  runs.  The  chosen  variable  mesh  with  mxn  =  16x32  =  512  nodes 
(Table  6.19)  over  the  half-length  and  half-circumference  is  based  on 
results  obtained  from  two  preliminary  linear  runs  in  which  the  mesh 
spacing  was  coarse  in  the  axial  direction(x)  and  fine  in  the  hoop 
direction(y)  for  the  first  run,  and  vice  versa  for  the  second. 
Approximately  five  elements  were  then  used  to  span  the  first  half 
wavelength  of  the  change  in  axial  curvature  in  the  hoop  and 
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axial  directions  obtained  respectively  from  the  first  and  second 
runs.+  The  size  of  consecutive  axial  and  hoop  intervals  was  then 
selected  to  increase  geometrically  by  20  and  10%,  respectively,  as 
may  be  deduced  from  Table  6.19.  Linear  results  appropriate  to  the 
finer  16x32  mesh  provided  about  the  same  wavelengths  as  those 
obtained  from  the  preliminary  runs,  and  were  found  to  agree  well  with 
available  series  solutions,  as  is  discussed  below.  Therefore,  the 
16x32  mesh  was  also  thought  to  be  adequate  to  describe  the  nonlinear 
behavior  of  the  cylinder  because  the  inward  deformation  pattern  under 
the  load  spreads  out  with  increasing  load,  and  this  results  in  an 
increasingly  greater  number  of  nodes  spanning  the  pattern. 

B.  Linear  Behavior  &  Comparison  with  Previous  Solutions 

Figures  6.36  and  6.37  show  the  hoop  variation  of  the  radial 
displacement  w  at  midlength  (x  =  0)  and  the  hoop  bending  moment  M 
at  the  axial  position  x  =  3.5" «( 1/140) (L/2) .  The  rapid  decay  of  the 
oscillatory  hoop  variations  with  increasing  circumferential  angle 
<t>  =  y/R  is  also  typical  of  the  STA6SC-1  results  for  the  axial 
bending  moment  and  the  axial  and  hoop  membrane  forces  Nx  and 
Ny.  In  contrast,  the  axial  variations  of  w  and  My  decrease 
monotonical ly  to  zero  at  the  simply  supported  edge  (x=L/2),  as  is 
shown  in  Figures  6.38  and  6.39  for  the  variation  of  w  along  the 
loaded  generator  ( $  =  0°)  and  along  a  neighboring  generator 
($  =  1/2°).  For  poked  shells,  note  that  the  slow  decay  of  w  in  the 
axial  direction  is  also  evident  from  graphs  given  in  [43]  for  longer 
cylinders  (L/R  =  4),  and  from  results  for  long  doubly  curved  shells 
presented  in  [44]  where  it  is  mentioned  that  displacements  spread 
much  further  in  the  direction  of  the  smaller  curvature  than  in  the 
direction  of  the  larger  one.  Also  note  that  the  axial  attenuation 
length  for  w  is  very  large  for  the  closely  related  case  of  a  pinched 
cyl inder  (see  [45]) . 


+*x  is  the  response  variable  with  the  shortest  significant  wavelengths 
in  both  directions.  Therefore,  results  are  more  accurate  for  the  more  slowly 
varying  and  lower  order  kinematic  variable  w(x,y),  the  radial  displacement. 
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The  above  linear  results  are  in  qualitative  agreement  with  series 

solutions  given  by  Bijlaard  [43,46,47]  and  Mizoguchi  [48]  for  thicker 

and  longer  cylinders.  Quantitative  comparisons  are  afforded  by  the 

series  solution  presented  by  Kempner,  Sheng,  and  Pohle  [49]  for  a 

range  of  shell  geometries  that  spans  the  geometry  (  R/t  =  638,  L/R  = 

.892)  of  the  poked  cylinder  analyzed  here.  The  radial  loading 

considered  in  [49]  is  an  axial  line  load  distributed  uniformly  over 

just  5%  of  the  shell  length,  and  located  symmetrically  about  the 

midlength  of  the  cylinder.  A  cylinder  with  such  a  short  line  load 

provides  a  good  check  case  for  the  point  loaded  cylinder.  Tabular 

results  are  presented  in  [49]  for  w,  M  ,  M  ,  N  and  N  at  the 

x  y  x  y 

following  points:(x,y)  =  (0,0),  the  location  of  the  resultant  load  P; 
(L/4,0),  the  “quarterlength"  point  on  the  loaded  generator  midway 
between  P  and  one  end  of  the  cylinder;  and  (0,L/4),  the  corresponding 
quarterlength  point  on  the  circumference  through  P.  Results  at  the 
first  point  for  R/t  =  500  and  800  and  L/R  =  .2  and  1  were 
logarithmically  interpolated  to  provide  values  of  the  radial 
displacement  and  the  stress  resultants  appropriate  to  R/t  =  638  and 
L/R  =  .892.  STAGSC-1  values  of  the  stress  resultants  at  the 
quarterlength  points  were  determined  from  known  values  at  neighboring 
centroidal+  points  by  means  of  linear  interpolation  in  one 
direction,  and  extrapolation  to  the  symmetry  plane  in  the  other 
direction.  The  extrapolated  values  were  obtained  by  passing  a  fourth 
order  polynomial  with  even  terms  (because  of  symmetry  considerations) 
through  three  points  and  by  evaluating  the  polynomial  at  the  origin. 
The  series  solution  [49]  and  STAGSC-1  results  for  the  three  points 
under  consideration  are  compared  in  turn  in  the  three  horizontal 
sections  of  Table  6.20.  Stress  resultants  that  are  small  or  which 


+The  option  to  print  out  values  of  the  stress  resultants  (or  stresses 
or  strains)  at  the  Gaussian  integration  points  doesn't  work.  Thus,  only 
centroidal  values  are  printed  out. 


0920B-86B:2 
( S3034)  30 


108 


I 


have  very  rapid  spatial  gradients+  at  the  quarterlength  points  are 
omitted  from  the  comparison.  Examination  of  the  table  reveals  that 
the  STAGSC-1  results  differ  from  the  series  solution  predictions  by 
an  average  of  only  +5.4%,  with  minimum  and  maximum  deviations  of  +.4 
and  +11.1%.  Furthermore,  each  STAGSC-1  result  is  above  the 
corresponding  series  solution  prediction,  and  thus  properly  reflects 
the  more  severe  nature  of  the  concentrated  loading  condition.  In 
view  of  the  interpolation  and  extrapolation  procedures  employed  to 
arrive  at  a  common  basis  of  comparison,  and  because  slight  phase 
shifts  or  dispersion  in  shape  at  points  of  moderate  spatial  gradients 
may  very  well  lead  to  percent  differences  on  the  order  of  those  given 
in  the  table,  it  seems  that  the  agreement  between  the  STAGSC-1  and 
series  solutions  for  the  two  slightly  different  loading  conditions  is 
about  as  good  as  can  be  reasonably  expected.  This  favorable 
correlation  of  linear  results  instills  a  sense  of  confidence  in  the 
reliability  of  the  STAGSC-1  code  in  general,  and  in  the  adequacy  of 
the  particular  mesh  selected  for  the  poked  cylinder  analysis, 
certainly  at  least  with  regard  to  the  prediction  of  linear  behavior. 

C.  Nonlinear  Behavior  and  Collapse 
Load-Displacement  Behavior 

The  continuous  load-displacement  curve  in  Figure  6.40  represents  the 
"best"  solution  obtained  here  by  refining  values  of  the  strategy 
parameters  in  the  DC  analysis.  Regions  I  and  II  in  the  figure 
represent  less  accurate  solutions  that  are  discussed  later.  The 
initial  linear  portion  of  the  curve  is  seen  to  extend  over  only  about 
10%  of  the  load  range  to  collapse  and  over  about  4%  of  the 
corresponding  displacement  range.  The  curve  then  bends  over  rapidly 


Comparisons  at  points  of  rapid  spatial  gradients  can  result  in  deceivingly 
large  percent  differences,  as  is  exemplified  by  consideration  of  two 
identical  shapes  that  are  slightly  out-of-phase  (shifted  relative  to  one 
another  in  the  coordinate  direction). 
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for  a  short  ranqe  of  load  levels  during  which  the  initially 
compressive  axial  and  hoop  membrane  stresses  near  the  load  become 
tensile  (see  Fiqure  6.41),  due  to  local  geometry  changes,  and  remain 
tensile  thereafter.  Figure  6.40  shows  that  collapse  occurs  at  the 
limit  point  of  the  curve  at  the  following  values: 

Pc  =  773,000  lb 
wc  =  65.6" 

or 


w  /t  =  37.5 

Somewhat  suprisingly,  the  post-collapse  behavior  could  not  be 
determined  from  the  DC  analysis.  Strategy  parameters  were  then 
varied  in  a  number  of  unsuccessful  runs,  each  one  being  terminated  at 
the  same  limit  point  with  a  message  that  the  stiffness  matrix  ceased 
to  be  positive  definite.  It  is  possible  that  the  load-displacement 
curve  is  very  steep  (stiffening  behavior)  immediately  after  the  limit 
point,  or  that  it  has  a  sharp  maximum  there,  or  both.  Another 
possible  explanation  is  given  after  considerat ion  of  spatial 
variation  plots. 


oat i a  1  Variations 


Hoop  and  axial  variations  of  w  and  at  collapse  are  displayed  in 
Figures  6.42  -  6.45  and  are  to  be  compared  with  the  corresponding 
spatial  variations  for  the  case  of  linear  behavior  at  low  load  levels 
(Figures  6.36  -  6.39).  Comparison  of  Figures  6.36  and  6.42  for  the 
hoop  variation  of  w  at  midlength  shows  that  the  inward  deformation 
pattern  under  the  load  spreads  out  with  increasing  load  (as  is  easily 
demonstrated  by  poking  at  a  thin-walled  beer  can)  and  that  the  shapes 
are  similar  otherwise.  In  contrast.  Figures  6.43  and  6.37  reveal 
that  the  hoop  variation  of  at  collapse  is  markedly  different 
from  that  at  low  load  levels.  The  jumpy  behavior  exhibited  in  Figure 
6.43  near  $  =  20°  is  due  to  the  development  of  a  sharp  hoop 
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curvature,  which  is  suggested  by  the  displacement  results  of  Figure 
6.42,  and  which  is  shown  directly  in  the  plot  of  the  change  of 
curvature,  given  in  Figure  6.46.  It  is  possible  that  the 
failure  to  trace  out  the  post-collapse  portion  of  the  P-wQ  curve 
may  be  due  to  ill-conditioning  associated  with  this  region  of  sharp 
curvature.  It  is  interesting  to  observe  from  Figure  6.46  that  the 
sharp  curvature  near  $  =  20°  is  preceded  by  a  spatial  region  in 
which  ity  «  -  l/R  so  that  the  total  curvature  +  1/R)  is 
approximately  zero,  or  in  other  words,  the  inward  displacement  U  < 
20°,  see  Figure  6.42)  is  circumferentially  flat  at  collapse  except, 
of  course,  in  the  immediate  neightborhood  of  the  load.  Figure  6.44 
shows  that  w  is  flat  also  axially  for  x>L/20.  Interestingly, 
comparison  of  the  load-deformation  plots  of  Figures  6.47  and  6.40  for 
iCy  near  $  =  20°  and  for  wq  reveals  that  starts  to  grow  rapidly 
at  a  load  level  (P/4  Z  150,000  pounds)  which  corresponds  to  the 
beginning  of  the  jog  (inflection  region)  in  the  P-wq  plot  of  Figure 
6.40. 

Scaled  Down  Model 


A  scaled  down  model  of  the  poked  cylinder  was  observed  to  deform  into 
a  spreading  diamond  shaped  pattern  with  rounded  "corners"  at 
midlength  (see  Figure  6.48).  The  circumferential  corners  became 
increasingly  sharper  with  increasing  load  until  a  subtle  and 
noncatastrophic  collapse  occurred  at  one  and  occasionally  both  of  the 
corners.  The  collapse  mechanism  appears  to  be  an  asymmetric  (about 
midlength)  local  snapping,  which  was  sometimes  accompanied  by  a  loud 
ping.  Such  collapse  behavior  was  noticed  at  most  but  not  all 
successively  poked  points  of  the  model.  Regardless  of  whether 
collapse  was  discernible,  stiffening  behavior  was  observed  to  take 
place  at  all  poked  points  after  the  circumferential  corners  became- 
sufficiently  sharp.  It  is  interesting  to  mention  that  these  corners 
appear  to  correspond  to  the  above  discussed  regions  of  sharp 
curvature  near  $  =  20°  at  midlength  that  are  predicted  by  the 
STAGSC-1  analysis. 


Ill 
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Mesh  Considerations  in  Retrospect 


As  mentioned  earlier,  the  chosen  mesh  was  thought  to  be  satisfactory 
for  the  nonlinear  analysis  since  it  provided  essentially  converged 
linear  results  that  are  in  good  agreement  with  available  analytical 
solutions,  and  because  the  number  of  nodes  spanning  the  spreading 
inward  deformation  pattern  under  the  load  increases  with  increasing 
load.  Evidently,  this  mesh  (or  any  finer  mesh)  is  inadequate  to 
describe  accurately  the  corner  near  $  =  20°  at  midlength.  However, 
it  is  important  to  point  out  that  the  solution  for  the  radial 
displacement  w(x,y)  may  still  be  accurate,  since  displacements 
converge  much  more  rapidly  than  higher  order  derivative  variables 
(My  or  Ky).  If  this  is  the  case,  and  since  the  P-wq  curve  of 
Figure  6.40  has  a  limit  point  (part  c  of  figure),  it  follows  that  the 


above  value  of  the  collapse 


load  at  this  point  would  also  be 


accurate. 


D.  Consideration  of  Strategy  Parameters 


It  has  already  been  shown  that  the  DC  solution  for  the  Venetian  blind 
drifts  away  from  the  more  accurate  LC  solution  in  a  local  region  of 
the  P-wo  curve.  Such  drift  was  also  observed  in  two  instances  in 
the  DC  solution  for  the  poked  cylinder,  as  shown  by  regions  I  and  II 
in  Figure  6.40.  Recall  that  the  continuous  curve  in  the  figure 
represents  the  "best"  solution  obtained  here  by  refining  values  of 
the  strategy  parameters  governing  the  convergence  criterion  (DELX) 
and  the  frequency  of  refactoring  (NEWT)  so  as  to  eliminate  drift. 
Specifically,  drift  in  region  I  was  eliminated  by  tightening  DELX  to 

-4  o 

10  from  its  default  value  of  10"  ,  and  drift  in  region  II  was 
smoothed  out  by  forcing  refactoring  at  every  step.  It  is  interesting 
to  note  that  the  drift  in  region  II  (see  Fig.  6.40b)  might  be 
interpreted  as  signifying  collapse,  which  would  be  spurious. 


In  view  of  the  improvement  brought  about  by  the  NEWT=-1  solution  in 
region  II,  an  additional  run  was  made  to  determine  if  the  drift  in 
region  I  could  similarly  be  avoided  by  refactoring  at  every  step, 
instead  of  by  tightening  the  condition  for  convergence.  That  this 
happens  is  evident  from  inspection  of  Figure  6.49  for  the  following 
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three  solutions,  which  are  shown  are  shown  to  an  expanded  scale:  (1) 
default  (DELX  =  10’^,  NEWT  t  -1),  (2)  tighter  convergence 
(10  \  /  -1)  and  (3)  successive  refactoring  (10’^,  -1).  It  is 
significant  to  observe  that  the  solution  based  on  the  more  stringent 
convergence  criterion  costs  2.25  times  as  much  as  does  the  one  which 
forces  refactoring  at  every  step,  since  it  requires  more  steps  (see 
Table  6.21).  Furthermore,  it  actually  requires  one  more 
refactoring,*  as  is  also  shown  in  Table  6.21.  Indeed,  the  solution 
based  on  successive  refactoring  turns  out  to  be  26%  cheaper  than  the 
default  solution,  since  it  required  roughly  half  as  many  steps. 
Similarly,  comparison  of  run  times  for  both  DELX=10-4  solutions  for 
region  II  of  Figure  6.40  shows  that  the  solution  with  successive 
refactoring  is  23%  cheaper.  It  appears  that  the  savings  due  to  the 
bigger  load  steps  allowed  by  successive  refactoring  more  than  offsets 
the  cost  of  the  additional  refactorings  for  this  size  problem. 

The  foregoing  results  now  make  it  possible  to  recommend  a  simple 
computational  procedure  for  the  economical  DC  analysis  of  the 
nonlinear  softening  behavior  of  shells  that  can  be  modelled  as  small 
or  modest  size  problems  on  STAGSC-1. 

E.  Recommended  Computational  Procedure  for  Displacement-Controlled 
Analyses 

Step  1:  Start  the  analysis  with  the  default  value  of  DELX=10~^  for 
the  convergence  criterion,  unless  there  is  a  special  reason 
to  proceed  otherwise,  and  let  the  program  decide  when 
refactoring  should  take  place  (NEWT>0), 

♦The  refactoring  that  always  occurs  during  the  first  step  of  a  restart  run  is 
excluded  from  this  relative  comparison. 
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Step  2:  Require  refactoring  at  every  displacement  step  ( NEWT = - 1 ) 

when  convergence  difficulties  or  inaccurate  solutions  are 
f irst  encountered.  Set  NEWT = - 1  for  the  rest  of  the 
analysis  until  stiffening  behavior  occurs  in  the 
post-collapse  portion  of  the  load-displacement  curve. 

Step  3:  Reduce  DELX  each  time  subsequent  computational  problems 
arise. 

Note  that  this  computat ional  procedure  may  not  be  economical  for 
problems  with  large  handwidths  because  of  the  cost  of  refactoring. 

The  present  results  indicate  that  it  is  economical  for  the  largest 
problem  considered  here,  which  has  512  mesh  points  and  an  average 
semi-bandwidth  of  128.  Also  note  that  the  procedure  does  not  apply 
to  a  load-controlled  analysis  since  such  an  analysis  generally 
requires  smaller  steps  for  satisfaction  of  convergence  than  does  a 
displacement-control  led  analysis.  This  is  illustrated  by  the 
Venetian  blind  results  (Figure  6.34). 

6.4.4  CONCLUDING  REMARKS  ANO  RECOMMENDATIONS 

The  quadrilateral  flat  plate  element  411  with  2x2  integration  points  was  used 
in  analyses  of  the  geometrically  nonlinear  behavior  and  collapse  of  a  point 
loaded  Venetian  blind  and  of  pinched  and  poked  cylinders.  Displacement- 
controlled  analyses  were  performed  for  all  three  shells  along  with  a 
load-controlled  analysis  of  the  Venetian  blind.  Conclusions  and 
recommendations  stemming  from  these  analyses  are  as  follows: 

A.  The  load  displacement  curves  up  to  the  limit  point  were  obtained  in  a 
straightforward  way  for  the  three  problems.  Changes  in  user  strategy 
during  the  course  of  the  analysis  were  required  only  for  the  poked 
cylinder  due  to  convergence  difficulties  or  indications  that 
inaccurate  solutions  had  been  accepted.  In  this  respect  it  is  worth 
noting  that  the  user  has  some  control  over  the  sophisticated  but 
flexible  strategy  in  STAGSC-l,and  this  allows  for  a  more  efficient 
use  of  the  program. 


■« 

,  *  . 
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B.  The  STAGSC-1  results  agree  well  with  linear  series  solutions  and  with 
a  limited  number  of  finite  difference  nonlinear  solutions. 

C.  A  displacement-controlled  (DC)  analysis  is  considerably  cheaper  than 
a  load-controlled  (LC)  analysis  of  the  same  problem,  but  may  require 
more  decisions  pertaining  to  the  selection  of  the  strategy  parameters 
governing  the  convergence  criterion  (DELX)  and  the  frequency  of 
refactoring  (NEWT).  Specifically,  the  DC  analysis  of  the  Venetian 
blind  was  6  times  cheaper  than  the  LC  analysis  and  required  8  times 
fewer  steps  to  reach  the  same  limit  point. 

D.  Because  satisfaction  of  the  same  convergence  criterion  is  so  much 
easier,  it  appears  that  too  large  a  step  is  sometimes  accepted  in  a 
DC  analysis.  This  may  cause  the  load-displacement  curve  to  drift 
away  from  the  "true"  solution  appropriate  to  the  given  mesh.  For  the 
Venetian  blind,  the  DC  curve  drifted  slightly  away  from  the  more 
accurate  LC  curve  in  a  local  region  and  then  properly  returned  to  it 
when  the  program  decided  to  refactor  at  a  later  step.  Drift  in  the 
poked  cylinder  results  required  a  restart  from  a  previous  accurate 
solution  with  different  values  selected  for  the  strategy  parameters. 
Thus,  drift  is  either  inconsequential  or  easily  corrected,  and  is  a 
small  price  to  pay  for  the  relatively  greater  economy  of  a  DC 
analysis. 

E.  For  modest  size  problems,  drift  can  be  eliminated  at  a  lower  cost  by 
refactoring  at  each  displacement  step  instead  of  sharpening  the 
convergence  criterion. 

F.  Based  on  these  results,  a  simple  computational  procedure  is 
recommended  for  the  DC  analysis  of  the  nonlinear  softening  behavior 
of  shells  that  are  run  as  small  or  modest  size  problems  on  STAGSC-1. 

G.  A  OC  load-displacement  curve  is  smoother  if  solutions  only  at 
refactored  steps  are  plotted. 


09208-868:2 
( S3034)  37 


H.  STAGSC-1  has  the  advantage  of  automatically  reducing  the  load  or 
displacement  step  when  convergence  difficulties  are  encountered,  such 
as  when  the  limit  point  is  approached.  Thus,  in  a  LC  analysis  it  is 
easier  for  STAGSC-1  to  home-in  on  the  collapse  load  to  any  degree  of 
accuracy  than  it  is  for  some  other  programs  in  which  input  values  of 
tiu*  load  step  cannot  ho  reduced  internally. 

I.  Since  knowledge  of  the  post-collapse  behavior  is  important, 
considerat ion  should  be  given  to  the  implementation  of  special 
techniques  [50]  that  allow  a  limit  point  to  be  passed  in  a  LC 
analysis. 

J.  The  scope  of  the  post-processor  (STAPL)  should  be  expanded  to  include 
variable-variable  plots  or  "history"  plots  to  give  plots  of  load  vs. 
displacement,  for  example,  and  to  include  "snapshots"  of  the  spatial 
variation  of  the  response  variable  along  a  selected  grid  line  at  a 
specified  load  level  (or  time  point). 

K.  Finally,  it  seems  appropriate  to  close  by  remarking  that  experience 
gained  here  engenders  a  sense  of  confidence  in  the  reliability  of  the 
elastic  nonlinear  capability  of  STAGSC-1. 

6.5  PROGRAM  EFFICIENCY 

The  evaluation  of  the  computational  efficiency  of  STAGSC-1  has  been  based  on 
execution  statistics  compiled  during  the  course  of  this  study.  It  is 
necessarily  subjective  because  all  the  work  has  been  performed  on  the 
Westinghouse  Power  Systems  Computer  Centre  COC  7600  installation.  This  system 
has  a  relative. y  small  central  core  memory  (SCM)  (151000g,  60  bit  words) 
which  is  augmented  by  600000g  words  of  large  core  memory  (LCM)  and  65-jq 
million  words  of  high  speed  disk  storage.  The  STAGSC-1  version  implemented  on 
this  system  is  not  configured  to  use  LCM  and  any  data  which  is  not  required  in 
core  is  buffered  to  and  from  disk. 

The  measure  of  cost  on  the  CDC  7600  is  the  computer  resource  unit  (CRU)  which 
is  determined  by  means  of  an  algorithm  combining  the  usage  of  small  and  large 
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core  memory,  mass  storage,  input/output  transfers,  number  of  disk  accesses  and 
central  processor  time.  The  dayfile  output  gives  a  breakdown  of  the  total  CRU 
usage  into  its  individual  components  thus  giving  insight  into  relative  costs 
of  computation  and  data  handling.  It  must  of  course,  be  emphasized  that  this 
is  subjective  because  of  the  weighting  assigned  by  the  charging  algorithm  to 
these  various  components.  Nevertheless,  such  statistics  do  shed  light  on  the 
performance  of  STA6SC-1  because  it  becomes  apparent  where  improvements  could 
be  made.  The  most  valuable  data  are  direct  comparisons  between  solutions 
obtained  by  STAGSC-1  and  other  programs  to  the  same  problem.  However,  such  an 
in-depth  investigation  would  have  exceeded  the  resources  available  for  this 
project  and  therefore  only  one  such  comparison  was  made. 

The  data  presented  cover  linear  and  nonlinear  static  analyses,  small  and 
larger  size  problems,  linear  and  nonlinear  dynamic  transient  and  nonlinear 
collapse  analyses.  Table  6.22  provides  details  of  the  execution  statistics 
for  fifteen  separate  analyses,  including  one  problem  using  the  MARC  program. 
The  central  processor  time  (CPU)  and  the  total  resources  used  (CRU)  are  given 
for  both  the  pre-processor  (STAGS1,  MARCPRE)  and  the  execution  phases  (STAGS2, 
MARCSTR ) .  For  all  but  the  smallest  of  problems,  the  resources  used  by  STAGS1 
are  a  very  small  proportion  of  the  total  and  will  not  be  discussed  further. 

The  final  column  of  Table  6.22  is  an  index  which  is  the  ratio  of  CRUs  to  CP 
hours  and  is  a  measure  of  the  execution  efficiency  (or,  more  realistically, 
inefficiency)  of  the  program.  Qualitatively,  it  is  the  ratio  of  the  total 
resources  used  to  the  resources  used  solely  in  computation.  The  average  value 
is  14.9  with  a  maximum  spread  of  -5.5  to  +6.0.  This  compares  well  with  the 
value  of  14.5  obtained  for  the  MARC  problem.  Table  6.24  presents  a  comparison 
between  several  structural  analysis  programs  in  use  on  the  Westinghouse 
system.  This  shows  that  STAGS  uses  the  most  total  resources  per  central 
processor  hour.* 

Returning  to  the  data  in  Table  6.22,  the  one  point  of  direct  comparison 
between  STAGSC-1  and  MARC  shows  that  although,  on  a  relative  basis,  MARC 


*An  interesting  feature  is  that  the  performance  of  two  different  versions  of 
ANSYS  is  apparently  very  different  and  is  a  consequence  of  the  later  version 
(Revision  3)  being  tailored  to  run  efficiently  in  the  interactive  mode. 
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performed  more  efficiently  than  STAGSC-1,  MARC  used  2.8  times  the  resources  to 
solve  the  same  problem.  Ultimately,  it  is  on  such  a  basis  that  a  program  must 
finally  be  judged  {other  things,  such  as  accuracy,  being  equal).  Table  6.23 
gives  details  of  the  way  in  which  the  computer  resources  were  used  in  a  number 
of  STAGSC-1  runs.  In  all  but  two  cases,  the  number  of  disk  accesses  was  the 
largest  contributor  to  resource  usage.  Although  this  is  a  direct  function  of 
the  particular  charging  algorithm  employed,  it  is  nevertheless  indicative  of  a 
potential  area  for  program  improvement.  On  the  CDC  7600,  this  would  be 
reduced  by  the  use  of  large  core  memory  (LCM). 

To  sum  up,  it  appears  that  STAGSC-1  uses  a  significantly  larger  amount  of 
total  resources  for  a  given  computational  effort  than  other  structural 
analysis  programs  on  the  Westinghouse  system.  Much  of  this  resource  usage  is 
due  to  accessing  disk  storage.  On  the  other  hand,  an  absolute  comparison 
between  STAGSC-1  and  MARC  suggests  that  STAGSC-1  is  much  more  efficient  in  its 
total  resource  usage  for  a  given  problem.  This  indication  is  supported  by 
user  experience  in  a  qualitative,  but  undocumented  sense. 
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TABLE  6.1 

LINEAR  CONVERGENCE  OF  ELEMENTS  410  and  420  FOR  THE 
CLAMPEO,  PRESSURE  LOADED  SQUARE  PLATE 


QUAF  410  QUARC  420 

No.  of 


Elements 

WA  (in) 

%  Error 

WA  (in) 

%  Error 

4 

.25103 

-8.78 

.24292 

-11.73 

9 

.26461 

-3.85 

.27104 

-1.51 

16 

.26969 

-2.00 

.27487 

-0.12 

26 

.27208 

-1.13 

.27579 

+0.21 

TABLE  6.2 

NONLINEAR  CONVERGENCE  OF  ELEMENTS  410,  411,  420,  422  FOR  THE 
CLAMPEO,  PRESSURE  LOADED  SQUARE  PLATE 

%  Error 

No.  of 


Elements 

QUAF  410 

QUAF  411 

QUARC  420 

QUARC  422 

4 

-13.6 

-5.51 

-12.60 

-6.75 

9 

-7.16 

-4.63 

-4.88 

-2.46 

16 

-4.75 

-3.38 

-3.32 

-1.67 

25 

-3.5? 

-2.69 

-2.54 

-1.45 

64 

-2.19 

-  - 
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TABLE  6.3 

NONLINEAR  FLAT  PLATE  ANALYSIS:  EXECUTION  STATISTICS 


Run 

Element 

Type 

No.  of 
Elements 

No.  of 
Steps 

Total 

CPU  secs 

CPU  secs 
/step 

Comments 

10 

QUAF 

410 

4 

19 

1.977 

0. 1041 

Step  halved  once. 

11 

QUAF 

410 

9 

19 

3.500 

0. 1842 

Step  halved  once. 
1  Refactoring. 

9 

QUAF 

410 

16 

19 

6.275 

0.3303 

Step  halved  once. 

1  Refactoring. 

12 

QUAF 

410 

2b 

1 1 

6.  106 

0.5551 

1  Refactoring. 

13 

QUAF 

411 

4 

11 

2.018 

0. 1835 

1  Refactoring. 

14 

QUAF 

411 

9 

11 

4.308 

0.3916 

1  Refactoring. 

15 

QUAF 

411 

16 

11 

6.974 

0.6340 

1  Refactoring. 

16 

QUAF 

411 

25 

11 

10.171 

0.9246 

1  Refactoring. 

17 

QUARC 

420 

4 

18 

3.749 

0.2083 

Step  halved  once. 

18 

QUARC 

420 

9 

19 

11.389 

0.5994 

Step  halved  once. 

19 

QUARC 

420 

16 

19 

17.885 

0.9413 

Step  halved  once. 

20 

QUARC 

420 

25 

11 

19.375 

1.7614 

1  Refactoring. 

21 

QUARC 

422 

4 

19 

5.056 

0.2661 

Step  halved  once. 
1  Refactoring. 

22 

QUARC 

422 

9 

11 

8.055 

0.7323 

1  Refactoring. 

23 

QUARC 

422 

16 

11 

14. 142 

1.2856 

1  Refactoring. 
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TABLE  6.3  (Continued) 


Run 

Element 

Type 

No.  of 
Elements 

No.  of 
Steps 

Total 

CPU  secs 

CPU  secs 
/step 

Comments 

24 

QUARC 

4  22 

25 

11 

21.934 

1.9940 

1  Refactoring. 

25 

QUAF 

64 

11 

13.456 

1.2232 

1  Refactoring. 

410 


0920B-64B: 2 
(S3034) 


121 


TABLE  6.4 

QUAF  410-STIFFNESS  EIGENMOOES 


Mode 

Eigenvalue 

Type 

Description 

1 

0.0 

Rigid  Body 

2 

0.0 

Rigid  Body 

All  of  the  "rigid  body"  modes  are  linear 
combinations  of  the  three  translational  and 

3 

0.0 

Rigid  Body 

the  three  rotational  modes.  However,  the 
rotations  about  the  normal  (Rw)  are  not 

4 

0.0 

Rigid  Body 

consistent  in  the  rigid  body  sense  with  the 
other  five  d.o.f.  There  appears  to  be  some 

5 

0.0 

Rigid  Body 

kinematic  freedom  in  the  shape  functions 
which  gives  rise  to  the  extra  rigid  body 

6 

0.0 

Rigid  Body 

freedom. 

7 

0.0 

Rigid  Body 

8 

1598 

Membrane 

Antisymmetric  biaxial  stretching  with 
quadratic  variation  normal  to  each  side. 

9 

1709 

Twist 

Pure  linear  twisting. 

10 

1923 

Bending 

Anticlastic  bending. 

11 

2111 

Bending 

Antisymmetric  bending  with  quadratic  and 
cubic  variation  of  normal  displacement. 

12 

2111 

Bending 

Same  as  Mode  11  with  sides  interchanged. 

13 

3571 

Bending 

Doubly  symmetric  spherical  bending. 

14 

41296 

Bending  + 
Twist 

Cubic  displacement  along  each  edge.  Diagonal 
bending  at  two  corners,  diagonal  twist  at 
the  other  two. 

15 

41296 

Bending  + 
Twist 

Same  as  Mode  14  with  corners  interchanged. 

16 

42260 

Membrane 

Quadratic  displacements  along  two  opposite 
sides  and  cubic  along  the  other  two. 

17 

42260 

Membrane 

Same  as  Mode  16  with  sides  interchanged. 

18 

2338060 

Membrane 

"Shearing"  mode  with  cubic  displacements 
along  all  four  sides. 

19 

2370460 

Membrane 

Biaxial  extension/compression  with  cubic 
displacement  along  each  side. 
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TABLE  6.4  {Continued) 


Mode 

Eigen  Value 

Type 

Description 

20 

2370460 

Membrane 

Same  as  Mode  19  with  extension/compression 
sides  interchanged. 

21 

2435990 

Membrane 

Similar  to  Mode  8,  but  with  the  emphasis  on 
the  corner  displacements  (as  opposed  to 
normal  rotations). 

22 

4285710 

Membrane 

Pure  biaxial  compression  (no  corner 
rotations). 

23 

6870160 

Bending  + 
Twist 

Diagonal  twist  at  each  corner.  Cubic  normal 
displacements  along  each  side.  Straight 
diagonals. 

24 

61861600 

bending  + 
Twist 

Similar  to  Mode  23,  but  with  anticlastic 
bending  of  the  diagonals. 
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TABLE  6.5 

SINGLE  ELEMENT  LOAD  CASES--MOOAL  ENERGY  DISTRIBUTION 


TABLE  6.6 

FLEXURAL  FREQUENCIES  IN  XZ-PLANE 

Frequency  (Hz) 

Mode  % 


No. 

STAGSC-1 

Exact 

Difference 

Comments 

1 

88.3277 

88.8164 

-0.6 

2 

544.409 

556.603 

-2.2 

Motion  in  the  XZ-plane  only. 
All  axial  and  torsional 

3 

1496.60 

1558.20 

-4.0 

freedoms  suppressed. 

4 

2869.83 

3054.05 

-6.1 

20  active  degrees  of  freedom 

5 

4631.57 

5048.56 

-8.3 

6 

6739.66 

7541.67 

-10.6 

7 

9140.57 

10533.4 

-13.2 

8 

11746.7 

14023.8 

-16.2 

9 

14379.4 

18012.7 

-20.2 

10 

16594.3 

22500.4 

-26.2 
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TABLE  6.7 

CLOSELY  SPACED  AND  EQUAL  MODES 


Frequency  (Hz) 


Mode 

% 

Dispt.  Ratio 

No. 

STAGSC-1 

Exact 

Difference 

(XZ/XY) 

Comments 

1 

80.3059 

80.7422 

-.54 

1 

107 

2 

81.1081 

81 .5496 

-.54 

107 

1%  difference  in 

cross-section 

1 

dimensions . 

3 

495.253 

506.003 

-2.13 

107 

4 

500.173 

511.063 

-2.13 

107 

5 

1362.69 

1416.54 

-3.8 

1 

107 

6 

1376.09 

1430.71 

-3.8 

io7 

1 

80.3058 

80.7422 

-.54 

6.66 

2 

80.3058 

80.7422 

-.54 

1 

6.66 

3 

495.253 

506.003 

-2.13 

1.148 

Equal  dimensions 

4 

495.253 

506.003 

-2.13 

1 

1.148 

5 

1362.66 

1416.54 

-3.8 

1 

1.041 

6 

1362.66 

1416.54 

-3.8 

1.041 
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TABLE  6.8 

FLAT  PLATE  CANTILEVER  FREQUENCIES 


Frequency  (Hz)  %  Error 

Mode  Mode 


No. 

Type 

Exact 

STAGSC-1 

MARC 

Zienk. 

STAGSC-1 

MARC 

Zienk. 

1 

First  Bending 

846 

750 

845 

826 

-11.3 

-0.12 

-2.4 

2 

First  Twist 

3638 

2172 

3651 

3728 

-40.3 

0.36 

2.5 

3 

Second  Bending 

5266 

3784 

5280 

5157 

-28.1 

0.27 

-2.1 

4 

Second  Twist 

11870 

6335 

12100 

12055 

-46.6 

1.9 

1.6 

TABLE  6.9 

1/4  WAVE  MODEL  FREQUENCIES 


Fundamental 

Mesh 

Frequency  (Hz) 

%  Error 

4x4 

81.58 

-2.73 

5x5 

82.45 

-1.69 

127 
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TABLE  6.10 

1/8  CYLINDER  MODEL  FREQUENCIES 


Frequency  Element 


Elements 

Mode 

n 

Hz 

%  Error 

Aspect  Ratio 

Comments 

8  Circ. 

1 

8 

75.79 

-9.6 

2.55 

3  modes  only 

2  Axial 

4 

10 

88.37 

-13.6 

requested 

5 

r> 

99.77 

-4.6 

16  Circ. 

1 

8 

74.81 

-10.8 

5.09  Circumferential 

2  Axial 

4 

10 

84.83 

-17.1 

refinement 

5 

6 

99.80 

-4.6 

7 

12 

108.15 

-22.3 

16  Circ. 

1 

8 

82.17 

-2.0 

2.55 

Axial  refinement 

4  Axial 

4 

10 

99.22 

-3.0 

5 

6 

103.63 

-0.9 

7 

12 

133.89 

-3.9 

16  Circ. 

1 

8 

84.15 

0.33 

1.27 

Axial  refinement 

8  Axial 

4 

10 

103.28 

0.92 

5 

6 

104.55 

-0.01 

7 

12 

141.69 

1.74 

TABLE  6.11 

1/2  CYLINDER  MODEL  FREQUENCIES 
Frequency  Hz 

STAGSC-1 


Mode  No. 

STAGSC-1 

Exact 

%  Error 

n 

Comments 

1 

82.17 

83.87 

-2.0 

8 

2 

82.39 

-- 

-- 

8 

Spurious  mode 

3 

86.35 

87.69 

-1.5 

7 

4 

86.63 

-- 

-- 

7 

Spurious  mode 
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TABLE  6.12 

1/2  CYLINDER  MODEL  FREQUENCIES 


Frequency  -  Hz 

Actual 


Mode  No. 

STAGSC-1 

Exact 

%  Error 

n 

1 

82.17 

83.87 

-2.03 

8 

2 

86.35 

87.69 

-1.53 

7 

3 

87.75 

89.92 

-2.41 

9 

4 

99.34 

102.3 

-2.89 

10 

5 

103.63 

104.6 

-0.93 

6 

TABLE  6.13 

CRITICAL  TIME  STEP  ESTIMATES  -  EXPLICIT  INTEGRATION 


Method 

Critical  At  (secs) 

Comments 

(i) 

4.94  x  10-6 

Based  on 

highest  frequency 

(ii) 

(a)  2.03  x  10-6 

(b)  4.08  x  10-6 

Users'  Manual  formulae. 

(iii) 

2.36  x  10-6 

Speed  of 

sound. 

0920B-86B:2 
( S3034)  52 


129 


TABLE  6.14 

SUMMARY  OF  RING  TRANSIENT  ANALYSES 
Time  Step 


Run 

Type 

Operator 

at  (sec) 

Comments 

10 

Elastic- 

Plastic 

Park 

2  x  10-6 

Final  time  reached  =  620  x  10"6  sec. 

Time  limit. 

1  1 

1  1  astir. - 
P 1 asl  ic 

Park 

Automat ic 

Initial  at  =  ?  x  10-6  SPC>  Increased 
to  16  x  10-6  sec.  Solution  breaks  down. 

12 

Elastic- 

Plastic 

Park 

2  x  10-6 

Final  time  reached  =  1002  x  10'6  sec. 
TAPE22  saved  for  restart. 

13 

Elastic- 

Plastic 

Park 

2  x  10-6 

Attempted  restart  failed.  Message  indi¬ 
cated  plasticity  data  could  not  be  found. 

14 

Elastic 

Park 

2  x  10-6 

Small  number  of  iterations  required  per 
time  step  (2  or  1). 

15 

Elastic 

Park 

20  x  10-6 

at  reduced  to  10  x  10'6  sec.  at 

Step  2.  Small  time  limit  (8  secs). 

16 

Elastic 

Park 

10  x  10-6 

2000  x  10"6  secs,  reached.  Increased 
at  caused  periodic  refactoring. 

17 

Elastic 

Trapezoidal 

10  x  10-6 

Periodic  refactoring. 

18 

t  last  ic 

Trapezo idal 

2  x  10-6 

No  refactoring.  Average  of  ?  iterations 
per  time  step. 

19 

Elastic 

Gear  2nd 
Order 

10  x  10-6 

Periodic  refactoring. 

20 

Elastic 

Gear  2nd 
Order 

2  x  10-6 

No  refactoring.  1  or  2  iterations  per 
time  step. 

21 

Elastic 

Gear  3rd 
Order 

2  x  10-6 

A  lot  of  refactoring  after  240  x  10“6 
sec.  Determinant  changed  sign  at  336  x 
10-6  sec. 

22 

Elastic 

Explicit 

1  x  10-6 

Divergence  occured  after  14  steps. 

23 

Elastic 

Explicit 

1  x  10-7 

Time  limit  after  265  steps. 

24 

Elastic 

Expl icit 

5  x  lO-7 

Divergence  occured  after  17  steps. 

26 

1  lastic 

txpl  ic  it 

1  x  lO'7 

Solution  obtained  out  to  1710  steps. 

Time  limit. 

130 
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TABLE  6.15 

NONLINEAR  TRANSIENT  EXECUTION  STATISTICS 
No.  of  CPU  CPU  secs/ 


Run 

Type 

Operator 

Steps 

secs. 

Step 

CRU 

CRU/Step 

Comments 

12 

Elastic- 

Plastic 

Park 

501 

200.45 

0.400 

0.798 

.00159 

Approx.  1  iteration/ 
step.  Also  plastic 
sub-iterations. 

14 

Elastic 

Park 

1000 

139.61 

0.140 

0.587 

.00059 

2  iterations/step. 

No  refactoring. 

16 

Elastic 

Park 

200 

45.88 

0.229 

0.169 

.00085 

2  to  7  iterations/ 
step.  Periodic 
refactoring. 

17 

Elastic 

Trapezoidal 

200 

63.28 

0.316 

0.215 

.00108 

2  to  6  iterations/ 
step.  Frequent 
refactoring. 

18 

Elastic 

Trapezoidal 

1000 

141.53 

0.142 

0.593 

.00059 

2  iterations/step. 

No  refactoring. 

19 

Elastic 

Gear  2nd 
Order 

200 

38.84 

0.194 

0.148 

.00074 

2  to  7  iterations/ 
step.  Periodic 
refactoring. 

20 

Elastic 

Gear  2nd 
Order 

1000 

119.04 

0.119 

0.511 

.00051 

1  or  2  iterations/ 
step.  No  refactoring. 

21 

Elastic 

Gear  3rd 
Order 

167 

44.85 

0.269 

0.156 

.00093 

2  to  20  iterations/ 
step.  Much 
refactoring. 

25 

Elastic 

Explicit 

1710 

316.54 

0.185 

1.271 

.00074 

No  iterations  or 
refactoring  necessary 

09208-868:2 
{ S3034)  54 


131 


Cable  6.16  Dimensions,  Material  Properties  and  Collapse  Loads 


SS  =  limply  support  (classical) 

F  =  free 

%  diff.  is  ralativa  to  STAGSA  collapse  load 


TABLE  6.17 


PARTICULARS  OF  FINITE  ELEMENT 
Average  Semi  Number  Number  of 


Problem 

Mesh 

DOF 

-Bandwidth 

of  Runs 

Steps 

to: 

i 

m 

P 

c 

P  . 
end 

P 

c 

Venetian  Blind 

un  if . 

,  10 

8 

6C3 

64 

DC:  1 

19 

30 

37 

LC :  2 

157 

157 

198 

P inched  Cylinder 

var . 

12 

8 

693 

61 

6 

45 

59 

248 

Poked  Cyl inder 

var . 

16 

32 

4185 

128 

11 

102 

102 

3717 

m,n  =  Number  of  mesh  points  in  axial  &  circumferential  directions,  respectively 

DOF  =  Number  of  active  degrees  of  freedom 

Pc  =  Collapse  load  at  limit  point  of  P-w  curve 

Pend  *  Denotes  end  of  analysis 

CPS  =  Computer  processing  seconds  on  the  Westinghouse  CDC-7600  computer 
CRU  =  Computer  Resource  Units 

$  Cost  is  based  on  in-house  rate  of  $90/CRU 

_3 

DELX  =  Error  tolerance  for  convergence.  Default  value  is  10 

(Note  that  the  default  value  of  the  relaxation  factor  WUND  was  used  in  all  anal 
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i 


lEFENT  ANALYSES 


|  CPS 


kefactcr  at 
Every  Step 


P  p  P 

end  c  end 


m 

65 

.170 

.297 

15 

27 

10-3 

No 

1 e 

198 

1 .021 

1 

.021 

92 

92 

10-3 

Nr 

18 

289 

.940 

1 

.098 

85 

99 

10-3 

Nc 

17 

3717 

21.367 

21 

.367 

1923 

1923 

10‘3,10'4,10'5 

Yes  for 

11  analyses) 


TABLE  6.18 

COMPARISON  OF  LINEAR  SOLUTIONS  FOR  VENETIAN  BLIND 
AND  PINCHED  CYLINDER  PROBLEMS 


Problem 

Mesh 

w  Et/P 

0 

(mxn) 

STAGSC-1 

Prev.  Sol. 

%  Diff.(a) 

Venetian  Blind 

10x8 

1091 

1156(b) 

-5.6 

Pinched  Cylinder 

12x8 

148 

164(c) 

-9.8 

23x15 

160 

-2.4 

(a)  %  difference 

is  relative 

to  the  previous 

solution 

(h)  Broqan  &  Almroth,  1971 
(c)  Lindherg,  et  al.,  1969 
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TA8LE  6.19 

VARIABLE  MESHES  USED  IN  PINCHED  AND  POKED  CYLINDER  ANALYSES 


Size  of  Intervals  in  Axial  &  Hoop  Directions 
Pinched  Cyl inder 


Row  (or  Col.) 


_ Pok ed  Cylinder 


Numbers 

Axial ( in. ) 

Hoop (deg)  Axial (in.) 

Hoop(deg) 

1-2 

.0124 

3.481 

7 

1 

2-3 

.0170 

5.469 

8.4 

1.1 

3-4 

.0232 

6.96 

10.08 

1.21 

4-5 

.0341 

10.94 

12.096 

1.331 

5-6 

.0449 

19.89 

14.5152 

1.4641 

6-7 

.0495 

21.38 

17.4182 

1 .6105 

7-8 

.0557 

21 .88 

20.9019 

1.7716 

8-9 

.0651 

25.0823 

1.9487 

9-10 

.0650 

30.0987 

2.1436 

10-11 

.0665 

36.1185 

2.3579 

11-12 

.0666 

43.3422 

2.5937 

12-13 

52.0106 

2.8531 

13-14 

62.4127 

3.1384 

14-15 

74.8952 

3.4523 

15-16 

83.6286 

3.7975 

16-17 

4.1172 

17-18 

4.5950 

18-19 

5.0545 

19-20 

5.5599 

20-21 

6.1159 

21-22 

6.7275 

22-23 

7.4002 

23-24 

8.1403 

24-25 

8.9543 

25-26 

9.8497 

26-27 

10.8347 

27-28 

11.9182 

28-29 

13.1100 

29-30 

14.4210 

30-31 

15.8631 

31-32 

15.5060 

Load  is  applied 

at  Row  1  and 

Column  1. 

The  size  of  consecutive  axial 

and  hoop  intervals  for 

the  poked  cylinder 

increase  geometrically  by  20% 

and  10%,  respectively, 

except  for 

the  last 

interval . 
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Table  6.20  Comparison  Of  Linear  Results  For  The  Poked  Cylinder 


Location  .  Variable  STAGSC-1  Kempner  et  al.  %0iff.<«) 


at  load,  x=y=*0 

w0  Et/P 

1092 

1045 

4.5 

x-L/4 

wEt/P 

369 

366 

.4 

y=0 

My/P 

.01542 

.01492 

3.3 

.00405 

.00395 

2.7 

NXR/P 

•15.04 

•13.72 

9.6 

x«0 

wEt/P 

346 

325 

6.5 

y=L/4 

nxr/p 

27.50 

24.75 

11.1 

_  I 

Average  =  5.4% 


(a)  %  Difference  is  Relative  to  Previous  Solution 
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TABLE  6.21 


COMPARISON 

OF  THREE  SOLUTIONS  FOR 

REGION  I 

Base  Run 

"Tighten 

"Refactored" 

(Default) 

Conv."  Run 

Run 

DELX  = 

10*3 

io-4 

IO’3 

Refactor  at 

No 

No 

Yes 

Every  Step? 

CRU  = 

2.582 

4.298 

1.908 

wQ(in. ) 

£  Ob) 

£  Ob) 

£  Ob) 

18.5 

69880+ 

69886+ 

19 

71592+ 

71523 

19.5 

73453 

73319* 

20 

75364 

20.5 

76888 

76626* 

21 

78449 

78360* 

21.5 

79771 

22 

82842* 

81538* 

22.5 

83104 

82933* 

23 

86392 

84742 

23.5 

86421* 

86113* 

24 

88421 

87593 

24.5 

90249 

89106 

88592* 

25 

90539* 

90460+ 

25.5 

92010 

91847* 

26 

92873 

93289* 

26.5 

“94727+ 

94826 

93625* 

27 

96537* 

27.5 

97309 

97359 

28 

98530 

28.5 

102340 

99774* 

98397* 

29 

101470 

29.5 

99486 

JSSSQ.Q 

30 

103870+ 

30.5 

105140* 

103020* 

Automatic  refactoring  during  first  step  of  restart  run 
♦Refactoring 
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TABLE  6.23 

STAGSC-1  EXECUTION  STATISTICS  -  2 


Percent  of  (otal  CPUs 


1/0  Uisk  Disk 


Run 

Description 

SCM 

I/O 

Buffers 

Storage 

Accesses 

CP 

CYLR00F 

05 

Gravity  loaded  cylindrical 
roof.  Linear  static. 

13.0 

17.4 

4.4 

0 

52.2 

13.0 

PLATETRAN 

06 

Linear  transient  analysis 
of  flatplate.  Trapezoidal 
integration 

13.2 

2.6 

2.6 

0 

73.7 

7.9 

MARC  DEMO 
PR0B401 A 

Same  as  above  using  MARC 
and  modal  superposition. 

16.1 

4.6 

1.1 

0 

67.8 

10.3 

RINGSTAT 

05 

Elastic-plastic,  static, 
large  disp.  Center 
loaded  ring. 

21.0 

22.6 

4.8 

0 

37.1 

14.5 

RINGTRAN 

16 

Elastic,  large  disp. 

Ring  with  initial 
velocity.  Park  method. 

15.9 

22.5 

3.3 

0.5 

47.3 

10.4 

RINGTRAN 

25 

Same  as  above  using 
explicit  method. 

15.1 

3.1 

2.1 

1.3 

68.4 

9.9 

WMSLOJA 

Pinch-loaded  cylinder. 
Nonlinear  collapse. 

16.0 

35.8 

2.8 

0.6 

34.3 

10.5 

WMSL07Y 

Point-loaded  cylinder 

Nonl inear  collapse. 

13.0 

42.9 

2.1 

1 .5 

33.8 

6.6 
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TABLE  6.24 

COMPARATIVE  PROGRAM  STATISTICS 


Program 

CRU 

CP  hr. 

STAGS 

16.6 

MARC 

11.5 

WECAN 

9. 1 

ANSYS 

Rev.  2 

5.8 

ANSYS 

Rev.  3 

12.0 

PLACRE 

10.6 

EDGE 
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,t  aurrunicu 


SHELL  WEIGHT  =  90.0  lb/ft2 

E  *  3.0  x  10®  Ib/in 
v  =0.0 


Figure  6. 1  Simply  Supported  Cylindrical  Roof 
-  Gravity  Load 


o  STAGSC  1  QUAF  410 
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Figure  6.3  Cylindrical  Shell  Roof  -  Element  Convergence  <  2 1 


DEFLECTION  w.  (in.) 


0.28 


0.27 


0.26 


0.25 


0.24 


0.23 


0.22 
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NUMBER  OF  ELEMENTS 


Figure  6.6  Flat  Plate  with  Clamped  Edges  Under  Uniform  Pressure  Loading 
Linear  Element  Convergence 
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PRESSURE  (Ib/iir) 


NUMBER  OF  ELEMENTS 


Figure  6.8  Flat  Plate  with  Clamped  Edges  Under  Uniform  Pressure  Loading 
Nonlinear  Element  Convergence 
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PERCENT  ERROR  IN  NONLINEAR  SOLUTION 


Figure  6.9  Flat  Plate  with  Clamped  Edges  Under  Uniform  Pressure 
-  Percent  Error  in  Nonlinear  Solution 
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Figure  6. 13  3-D  Beam  Model  For  Modal  Analysis 
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Figure  6.20  1/2  Cylinder  Model  -  Mode  Shapes  (2l 
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Figure  6.24  Cantilever  Plate  Under  Transient  Pressure  Loading 
-  Tip  Displacement  History 
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Figure  6.27  Impulsively  Loaded  Ring 

-Nonlinear  Elastic  Plastic  Response 
—Deflection 


170 


CIRCUMFERENTIAL  STRAIN,  OUTER  SURFACE,  RING  CENTER 


Figure  6.29.  Impulsively  Loaded  Ring  •  Deformed  Shape  at  About  0.0008  Seconds 
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Figure  6.3i  Impulsively  Loaded  Ring  •  Nonlinear  Elastic  Response 
—  Park  Operator 
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Figure  6.32  Impulsively  Loaded  Ring  -  Nonlinear  Elastic  Response 
—  Trape/.oidal  Operator 
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Figure  6.33  Impulsively  Loaded  Ring  -  Nonlinear  Elastic  Response 
-  Gear's  Operator 
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Figure  6.34  Comparison  of  Load-Displacement  Results  for  the  Venetian  Blind 
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f  igure  6.40  Load-Displacement  Results  fur  Poked  Cylinder 
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c)  Rtfion  ill 


Figure  6.40  (Continued) 
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h igure  6.46.  Hoop  Variation  Along  \  =  3.5"ol'  Change  in  Hoop  Curvature  at  C  ollapse. 
Poked  Cylinder 


5407-96 


190 


(qi)  r.Ol  X 


Ky  (1/in.) 


Figure  6.47.  Load  vs.  Change  in  Hoop  Curvature  at  x  =  3.5".  <*>  -  22.95°.  Pok 
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7.0  CONCLUSIONS  AND  RECOMMENDATIONS 


In  the  Introduction,  it  was  pointed  out  that  the  criteria  for  evaluating  a 
structural  analysis  computer  program  should  reflect  both  the  nature  of  the 
program  and  the  class  of  users  most  likely  to  require  the  program.  Since 
STAGSC-1  is  a  nonlinear  shell  analysis  program  its  nature  is  quite  specialized 
and,  by  the  same  token,  the  majority  of  users  will  be  relatively 
sophisticated.  This  places  considerable  weight  on  the  technical  excellence  of 
the  program,  its  flexibility  in  use  and  for  post-processing  and  finally,  on 
its  documentation. 

In  order  to  obtain  a  systematic  evaluation  of  the  different  aspects  studied,  a 
rating  will  be  assigned  to  each  aspect  based  on  a  scale  of  0  to  4.  Table  7.1 
shows  the  interpretation  of  the  rating  method. 

For  the  purposes  of  rating,  five  major  aspects  have  been  identified  and  given 
a  relative  weight  as  follows: 


0 

Program  documentation 

20% 

0 

Analysis  capability 

20% 

0 

Input  and  output 

20% 

0 

Program  performance 

25% 

0 

User  support 

15% 

These  will  each  be  discussed  separately  in  the  next  sub-section,  but  at  this 
point  the  overall  conclusions  can  be  stated;  STAGSC-1  rates  as  about  mid-way 
between  "acceptable  and  "good"  (i.e.,  an  arithmetic  result  of  2.39). 

This  may  seem  to  be  a  somewhat  conservative  evaluation  considering  the  many 
unique  capabilities  of  the  program.  In  the  opinion  of  the  authors  of  this 
report,  however,  the  lack  of  documentation  of  the  program  itself  and  the 
shortcomings  in  the  area  of  post-processing  justify  this  conclusion.  If  these 
two  aspects  were  each  upgraded  to  "excellent",  then  on  the  present  system  of 
rating  the  overall  evaluation  would  than  place  STAuSC-1  in  the 
"good-to-excel lent"  range. 
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7.1  PROGRAM  DOCUMENTATION 


The  three  types  of  documentation  evaluated  were  the  users'  manual,  theoretical 
manual  and  programmers'  manual. 

A.  Users'  Manual  (Rating  3.5) 

Highly  rated  because  of  logical  organization  and  excellent  user 
guidance.  Negative  features  are  lack  of  any  index  or 
cross-referencing  and  insufficiently  clear  indication  of  currently 
inoperational  features. 

B.  Theoretical  Manual  (Rating  1.5) 

The  low  rating  is  due  to  incompleteness  and  substantial  obsolescence. 

C.  Programmers'  Manual  (Rating  0) 

The  programmers'  manual  does  not  exist.  Equal  weighting  is  given  to 
each  of  these  so  that  the  overall  rating  for  documentation  is  only 
1.67,  i.e.,  between  "maryinal"  and  "acceptable". 

7.2  ANALYSIS  CAPABILITY 

Eleven  areas  were  selected  for  the  purposes  of  rating.  Each  will  be  briefly 
discussed  and  rated. 

A.  Static  Linear  (Rating  3.5) 

STAGSC-1  is  not  specifically  intended  for  linear  analysis  but  such 
problems  are  nevertheless  handled  very  well  by  the  program. 
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B.  Static  Nonlinear  (Rating  4.0) 


The  highest  rating  must  be  given  to  this  capability  by  virtue  of  its 
speed  and  the  way  in  which  solution  strategy  is  handled 
semi -automat ical ly  by  the  program  while  leaving  the  user  significant 
control . 

C.  Vibration  Modes  (Rating  3.5) 

This  feature  is  also  highly  rated  because  of  the  use  of  a  very 
effective  eigenvalue  solver,  i.e.,  subspace  iteration.  This 
capability  is  current  "state-of-the-art." 

D.  Bifurcation  Buckling  (Rating  3.5) 

The  same  comments  apply  as  in  the  case  of  vibration  analysis. 

E.  Plasticity  (Rating  2.0) 

The  plasticity  capability  is  judged  as  "acceptable"  for  the  reason 
that  although  a  relatively  modern  hardening  model  (White-Bessel ing ) 
is  implemented,  current  practice  requires  that  a  variety  of  hardening 
models  be  available.  Therefore,  as  a  minimum,  isotropic  hardening 
should  also  be  implemented  (kinematic  hardening  is  available  as  a 
degenerate  form  of  White-Bessel ing) . 

F.  Nonlinear  Collapse  (Rating  3.8) 

This  capability  only  falls  short  of  the  "excellent’'  category  because 
such  analysis  will  often  require  the  use  of  plasticity  with  the 
consequent  limitations  as  described  above. 
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G.  Transient  Response  (Rating  3.8) 


The  time  integration  procedures  available  in  STAGSC-1  provide 
standard  explicit  (central  difference)  and  implicit  (trapezoidal) 
methods  plus  three  other  implicit  schemes  (Park  and  Gear's  2nd  and 
3rd  order).  There  is  therefore  considerable  scope  for  the  analyst  to 
select  a  procedure  best  suited  to  the  problem  in  hand. 

H.  User  Subroutines  (Rating  4.0) 

The  wide  variety  of  user-coded  functions  which  can  be  specially 
provided  for  any  given  function  is  well  suited  to  the  sophisticated 
user.  The  range  of  capabilities  available  in  this  form  is  believed 
to  be  exceeded  only  by  the  MARC  Program  [26]. 

I.  Element  Library  (Rating  1.5) 

The  heart  of  any  finite  element  program  is  the  library  of  elements. 

In  the  case  of  STAGSC-1,  the  library  consists  mostly  of  flat  plate 
elements  with  membrane  and  bending  capability.  There  are,  however, 
no  genuine  shell  elements  with  curvature  and  this  must  be  considered 
a  significant  disadvantage  in  a  shell  analysis  program. 

J.  Elastic  Material  Properties  (Rating  3.0) 

Both  isotropic  and  orthotropic  material  behavior  are  provided  for 
together  with  the  ability  to  orient  material  axes  in  any  given 
direction.  Multilayer  composite  materials  may  be  represented  by 
specifying  properties  for  each  layer  independently. 

K.  Temperature  Effects  (Rating  1.0) 

Temperatures  may  be  specified  as  a  function  of  surface  coordinates 
but  not  in  the  thickness  direction  (except  for  discrete  stiffners). 
Moreover,  material  properties  cannot  be  expressed  directly  as  a 
function  of  temperature. 
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The  overall  rating  for  analysis  capability  is  2.89  (i.e.,  "good")  with  the 
element  library  carrying  a  weighting  of  3  and  all  other  components  weighted  as 
1  or  less. 

7.3  INPUT  AND  OUTPUT 

A.  General  Input  (Rating  3.5) 

The  free  format  style  of  input  and  the  ability  to  intersperse  comment 
statements  make  for  ease  of  input  and  subsequent  ease  of  reading. 

The  volume  of  input  required  is  not  excessive,  even  for  larger 
problems . 

B.  Shell  Unit  Input  (Rating  3.5) 

The  ability  to  generate  both  regular  and  irregular  mesnes  in  standard 
shell  surface  geometries  is  an  excellent  feature.  Shell  connections 
are  also  easily  specified  if  the  shell  units  are  joined  at  compatible 
boundaries  (e.g.,  a  cylinder  to  a  torus  at  a  cross-sectional 
boundary).  More  general  intersections  are  not  determined  by  tne 
program  (e.g.,  a  cylinder/sphere  intersection)  and  the  intersecting 
boundaries  must  be  calculated  and  input  by  the  user.  This  also 
denies  the  use  of  the  shell  unit  library  and  the  intersecting  shells 
must  be  defined  by  user  subroutine  or  by  a  spline  fit  to  specified 
points. 

C.  Element  Unit  Input  (Rating  1.5) 

Input  for  the  element  unit  places  a  much  greater  burden  on  the  user 
and  the  generation  of  a  finite  element  mesh  requires  the  use  of  user 
subroutines;  otherwise,  input  must  be  provided  node  by  node  and 
element  by  element.  However,  the  use  of  the  element  unit  is  normally 
expected  to  be  confined  to  local  regions  of  the  shell  where  the 
geometry  and/or  the  mesh  may  be  very  irregular. 
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D.  User  Subroutines  (Rating  4.0) 


As  has  been  indicated,  most  of  the  input  data  may  be  specified  by 
user  subroutines  and  in  many  instances  they  provide  the  only  means 
(e.g.,  shell  geometries  not  included  in  the  library  or  shell  surface 
imperfections! .  The  capability  must  be  rated  as  excellent  on  account 
of  the  range  ot  functions  available. 

E.  Gc  /era!  Output  (Rating  2.0) 

There  is  an  acceptable  degree  of  selective  control  over  solution 
output,  at  least  in  terms  of  the  frequency  of  the  output.  According 
to  the  manual,  complete  sets  of  displacements  or  stresses,  stress 
resultants  etc.  can  be  printed  at  different  intervals;  this  feature 
is  currently  defective  in  that  all  the  various  quantities  must  be 
output  at  every  load  or  time  step.  A  frequency  specification  for 
selected  displacements  etc.  would  be  a  considerable  improvement. 

The  occurrence  of  yielding  is  noted  in  the  stress  and  strain  output 
but  only  the  total  strains  are  given.  There  should  also  be  an  option 
to  obtain  the  elastic  and  plastic  components  of  the  total  strain. 

F.  Post-Processing  (Rating  0.5) 

The  currently  available  post-processing  using  the  STAPL  program  is 
inadequate  for  many  analysis  tasks.  The  model  and  solution  data  are 
saved  for  restart  and  this  should  be  made  accessible  to  the  user  tor 
post-processing  in  any  desired  manner.  A  minimum  effort  in  this 
direction  would  be  for  the  developer  to  provide  a  complete 
description  of  the  model  geometry  and  solution  data  files. 

The  overall  rating  for  input  and  output  is  2.50  with  equal  weighting  given  to 

each  topic  discussed. 
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7.4  PERFORMANCE 


The  evaulation  of  performance  is  based  ori  the  specific  advanced  evaluation 
studies  described  in  Section  6.  Since  this  examines  only  selected  aspects  of 
the  performance  of  STAGSC-1  the  results  cannot  be  claimed  to  be  comprehensive 
but  are,  at  least,  a  reasonable  indication  of  general  performance. 

A.  Element  Convergence  (Rating  1.5) 

The  element  convergence  observed  in  the  tests  conducted  is  judged  to 
be  less  than  acceptable  in  comparison  with  “state-of-the-art"  curved 
shell  elements.  This  is  regarded  as  a  serious  shortcoming  not  only 
because  of  the  loss  in  accuracy  but  also  because  the  full  potential 
of  STAGSC-1  as  an  efficient  shell  analysis  program  cannot  be  properly 
exploited.  The  promised  implementation  of  the  Ahmad  isoparametric 
shell  element  will,  it  is  hoped,  greatly  improve  this  crucial  aspect 
of  performance. 

B.  Eigenvalue  Solutions  (Rating  3.5) 

The  eigenvalue  solver  performed  well  in  determining  modes  and 
frequencies  in  the  often  pathological  situations  of  closely  spaced  or 
repeated  eigenvalues.  It  was  also  seen  to  be  effective  in  computing 
a  mode  closest  to  a  specified  frequency.  Situations  where  the 
performance  was  less  effective  were  seen  to  be  a  consequence  of  an 
unsatisfactory  mesh  or  element  aspect  ratio  and  hence  not 
attributable  to  the  eigensolver. 

C.  Transient  Integration  (Rating  3.5) 

Performance  of  the  transient  integration  schemes  for  linear  and 
highly  nonlinear  problems  was  generally  good.  The  trapezoidal  scheme 
appeared  marginally  the  best  based  on  the  least  amount  of  frequency 
distortion  and  damping  for  the  nonlinear  transient.  The  explicit 
method  was  unable  to  complete  the  nonlinear  problem  even  with  an 
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order  of  magnitude  smaller  time  step  than  the  estimated  critical 
value.  From  the  aspect  of  economy,  the  explicit  method  did  not 
appear  advantageous  because  the  problem  size  was  too  small.  However, 
for  highly  nonlinear  problems,  the  implicit  schemes  are  of  greater 
utility  since  the  stiffness  matrix  can  be  reformulated. 

D.  Nonlinear  Collapse  (Rating  3.5) 

The  most  significant  conclusion  that  may  be  drawn  from  the 
performance  of  the  elastic  nonlinear  collapse  analyses  is  that 
STAGSC-1  has  a  clearly  demonstrated  capability  to  trace  out  a 
load-displacement  curve  for  a  shell  with  somewhat  subtle  collapse 
behavior.  In  the  case  of  the  point  loaded  cylinder,  the  predicted 
deformation  pattern  up  to  collapse  appears  to  be  in  qualitative 
agreement  with  that  observed  visually  in  a  small  scale  model  of  the 
cyl inder . 

STAGSC-1  performed  successfully  both  for  controlled  load  and 
controlled  displacement  incrementation,  largely  because  of  the 
adaptive  incrementation  based  on  the  convergence  of  the  solution. 

This  particular  feature  distinguishes  STAGSC-1  from  most  other  finite 
element  programs  with  nonlinear  capability.  Sufficient  control  was 
found  to  be  available  to  the  user  both  for  incrementation  strategy 
(ICUT  and  NEWT)  and  convergence  tolerance  (DELX).  In  addition  to 
effective  control  of  solution  accuracy,  these  strategic  parameters 
provide  the  means  to  obtain  the  most  economical  solution. 

E.  Program  Efficiency  (Rating  3.0) 

The  overall  efficiency  of  STAGSC-1  cannot  be  adequately  assessed  in 
relation  to  other  programs  based  on  the  present  study.  However, 
given  the  available  evidence  it  appears  that: 

(1)  STAGSC-1  is  a  factor  of  2  to  3  times  more  economical  in  solving 
a  small  linear  transient  problem  than  the  MARC  program. 
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(Z)  there  appears  to  be  scope  for  improvement  in  data  management  for 
the  COC  7600  system  version.  Based  on  the  admittedly  subjective 
and  imprecise  impressions  of  other  users,  5 IAG5L-1  is  a 
relatively  fast  running  program  in  its  present  form.  It  seems 
safe  to  conclude,  therefore,  that  improvements  in  the  element 
library  and  in  data  handling  could  bring  improvements  in  overall 
efficiency  which  would  make  it  a  first  choice  for  shell  analysis 

Giving  equal  weighting  to  each  of  the  performance  aspects  of  yields  an  overall 
rating  of  3.0. 

7.5  USER  SUPPORT  (RATING  1.5) 

The  topic  of  developer  support  for  users  of  STAG5C-1  has  not  been  discussed 
elsewhere,  but  being  a  significant  factor  in  the  use  of  a  large  scale  program 
it  deserves  some  consideration.  STAGSC-1  is  a  program  in  the  public  domain 
and  the  developers,  Lockheed  Palo  Alto  Research  Laboratory,  are  not  suppliers 
of  structural  software  in  a  general  commercial  sense.  Potential  users  can 
obtain  the  source  file  for  a  nominal  fee  and  are  then  responsible  for 
installation  on  their  own  systems.  For  an  additional  annual  fee,  updates  and 
consulting  services  are  available  from  Lockheed. 

Developer  expertise  is  available  for  solving  problems  and  applying  the  program 
correctly.  Program  malfunctions  can  be  reported  but  corrective  action  must  be 
taken  by  the  user  based  on  program  "fixes"  supplied  by  the  developer. 

Since  funding  for  STAGSC-1  development  is  principally  from  NASA,  the 
individual  user  has  little  or  no  control  over  future  development  items.  It  is 
clear  that  much  better  user  support  could  be  obtained  if  the  fees  were 
adequate  to  support  a  group  whose  sole  responsibility  was  consultation  and 
program  maintenance.  The  present  arrangement  fails  to  be  satisfactory  if  a 
user  requires  a  minor  modification  for  his  particular  purpose.  Generally 
speaking,  the  complexity  of  STAGSC-1  makes  user  written  modifications  to  the 
program  very  difficult  to  achieve  and  there  seem  to  be  no  proper  channels 
through  which  he  can  obtain  them  from  the  developers. 
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As  far  as  program  updates  are  concerned,  none  have  been  received  since  the 
evaluation  version  was  delivered  in  November  1979. 

During  the  present  evaluation,  the  staff  at  Lockheed  were  consulted  on 
numerous  occasions  and  excellent  cooperation  was  obtained. 

7.6  RECOMMENDATIONS 

A.  Documentation 

There  are  four  major  recommendations; 

(1)  An  up-to-date  theoretical  manual  with  adequate  coverage  of 
element  formulations  should  be  provided. 

(2)  A  programmer's  manual  is  required  which  describes  the  program 
flow  in  sufficient  detail  to  enable  the  skilled  user  to  identify 
the  source  of  any  program  difficulties  and  to  make  changes  with 
help  from  the  developers. 

(3)  The  users'  manual  should  be  provided  with  a  suitable  index. 

(4)  A  problem  demonstration  manual  should  be  written  specifically 
for  STAGSC-1 . 

8.  Capabilities 

(1)  A  good  curved  shell  element  (e.g.,  the  Ahmad  isoparametric)  is 
essential  for  the  full  potential  of  STAGSC-1  to  be  realized. 

The  simpler  elements  should  be  retained  (QUAF410,  411)  but  the 
value  of  the  higher  order  (QUARC)  series  would  then  be  less  and 
they  could  probably  be  discarded. 
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(2)  The  option  for  additional  plastic  hardening  models  should  be 
provided.  As  a  minimum,  isotropic  hardening  should  be 
implemented  and  consideration  should  be  given  to  simulating 
cyclic  hardening  (or  softening),  e.g.,  the  so-called  ORNl  10th 
Cycle  Hardening  Rule  [40]. 

(3)  Inelastic  behavior  should  be  extended  to  include  creep.  A 
simple  creep  law  (e.g.,  the  Norton  formula)  should  be  provided 
with  user  input  coefficients  and  the  option  for  time  hardening 
or  strain  hardening.  In  addition,  provision  should  be  made  for 
a  user  specified  creep  law  to  be  input.  Consideration  should 
also  be  given  to  the  eventual  inclusion  of  cyclic  creep  strain 
hardening  behavior,  e.g.,  the  so-called  ORNL  auxiliary  creep 
hardening  rules  [40]. 

(4)  Temperature  distributions  should  include  variations  through  the 
shell  wall  thickness.  In  addition,  material  properties  should 
be  allowed  to  vary  with  temperature. 

(5)  The  input  for  the  element  unit  could  be  improved.  Simple 
pattern  generation  for  nodal  coordinates  and  element 
connectivity,  in  addition  to  the  present  node-by-node  and  user 
subroutine  methods,  would  make  this  feature  much  more  usable. 

(6)  The  selectivity  of  stress  and  strain  printout  should  be  extended 
to  include  the  ability  to  print  selected  guantities  at  any 
desired  load  or  time  step  intervals.  Also,  the  option  should  be 
provided  to  obtain  strain  components  decomposed  into  elastic, 
plastic  and  thermal  contributions. 

(7)  A  user  post-processing  file  should  be  created  which  contains  all 
the  necessary  model  geometry  and  solution  data  to  permit 
solution  processing  according  to  the  needs  of  the  analyst.  The 
file  should  be  constructed  so  that  the  information  retrieval  is 
straightforward  and  it  should  be  fully  documented. 
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TABLE  7.1 

EVALUATION  RATINGS 


Rating 


Interpretation 


4 

3 

2 

1 

0 


Excellent 

Good 

Acceptable 

Marginal 

Poor 
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9.0  APPENDIX 


9.1  STAGSC-1  SUBROUTINE  LINKAGES 

The  overlay  structure,  as  described  In  Section  3,  Program  Architecture,  makes 
frequent  reference  to  the  various  subroutine  calls  made  In  each  overlay.  The 
diagrams  showing  the  subroutine  linking  within  each  overlay  are  contained  in 
this  section.  They  are  not  flow  diagrams  since  there  is  no  sequence  implied 
in  the  way  the  diagrams  are  constructed.  Some  of  the  overlays  are  too  large 
to  fit  all  the  subroutine  calls  within  a  single  diagram,  Fpp  example.  Overlay 
(1,0)  requires  three  diagrams  (Figures  9.4,  9.5,  and  9,6).  Where  there  are 
links  between  subroutines  in  different  figures,  the  connection  to  the  calling 
subroutine  is  shown  by  a  dashed  line. 

There  are  also  two  subroutines  which  make  a  considerable  number  of  calls  and 
are  themselves  basic  modules.  These  are  MACUP  in  STAGS1,  and  CVR1  and  VRDATA 
in  STAGS2.  The  details  of  the  calls  within  these  subroutines  are  given  only 
once  (Figure  9.11  for  MACUP  and  Figure  9.27  for  VRDATA). 

9.2  ELEMENT  STIFFNESS  EIGENSOLUTION 

The  concept  of  stiffness  matrix  elgenmodes  is  discussed  by  Gallagher  [29].  In 
order  to  obtain  the  distribution  of  strain  energy  with  respect  to  the 
eigenmodes  for  a  given  loading  case,  it  Is  only  necessary  to  decompose  the 
solution  vector  into  modal  components.  The  calculation  of  the  strain  energy 
is  then  straightforward  and  can  be  obtained  mode  by  mode  since  the  modes  are 
orthogonal.  The  necessary  equations  will  now  be  developed. 

For  an  element  stiffness  matrix  [k],  the  characteristic  polynomial  may  be 
obtained  by  expanding  the  determinant 

10]  -  x  [I]  I  *  0  9.1 

and  the  eigenvectors  by  solution  of  the  equation 
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03  (di }  =  x.  { } 


9.2 


where  {d^},  a.  are  the  i  th  eigenvector  and  eigenvalue. 

Because  of  the  orthogonality  of  the  eigenvectors,  we  may  write 

{d^1  [kj  {di }  =  \i  9.3 

provided  that 

{di }l  {d. }  =  1  9.4 

Thus,  the  eigenvalues  are  the  corresponding  generalized  stiffness  coefficients. 
A  matrix  of  the  eigenvectors  [rd]  may  be  defined  as 

[rd ]  =  [(d.j }  ....  { d ^ }  ....  {d^}]  9.5 

where  n  is  the  order  of  [k ]. 

Since  the  eigenvectors  are  linearly  independent,  any  solution  vector  {a}  may 
be  expressed  as  a  linear  combination  of  the  eigenvectors  as  follows: 

< A}  =  [rd]  {a}  9.6 

where  (a)  is  a  vector  of  modal  coefficients. 

Thus, 

la)  ■  [rd]_1  {4}  9.7 

which,  again  because  of  orthogonality,  may  be  written  equivalently  as 
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9.8 


U>  =  t^]1  (  a  } 

The  strain  energy  associated  with  the  displacement  vector  (a)  is 

U  =  \  {a}1  [k]  {*}  9.9 

Because  of  Equation  9.6,  this  may  also  be  written  as 

U  =  {a}1  [i^]1  [k]  [rd]  {a}  9.10 

Because  of  the  orthogonality  as  expressed  by  Equation  9.3,  this  may  be 
rewritten  as 

U  =  2  l®)*  A 
1  n  2 

=  •*  z  x-  9.12 

1  i  =  l  1  1 

where  A  is  the  diagonal  matrix  of  the  eigenvalues. 

The  program  ELMOD  was  therefore  written  to  perform  the  eigensolution  for  the 
element  stiffness  matrix  and  also  to  determine  the  corresponding  modal 
coefficients  for  a  solution  vector  (a).  The  strain  energy  distribution  is 
then  obtained  simply  from  Equation  9.12.  Table  9.1  gives  the  listing  for 
ELMOD. 
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Table  9.1  ELMOD  Listing 


PROGRAM  ELMOD ( INPUT , OUTPUT , TAPE5- INPUT , TAPE6=0UTPUT , TAPE10 , TAPE11 ) 
DIMENSION  A(1400) , ARE (50, 50) , EIV(50, 50) , EIG(50) ,D(50) , ALP (50) , 
1U(50) ,FV1(50) ,FV2(50) 

INTEGER  EIGEN, ATRANB 
EIGEN  =  3 
ATRANB  =  23 
C 

C  READS  IN  ORDER  OF  ELEMENT  MATRIX  (M) , FLAG  ITAP  TO  INDICATE  WHETHER 

C  STIFFNESS  DATA  AND  SOLUTION  VECTOR  ARE  TO  BE  READ  FROM  TAPE  OR 

C  CARDS . 

C 

C  ITAP  =  C  DATA  FROM  CARDS 

C  ITAP  =  1  DATA  FROM  TAPE 

C  NSKIP  IS  THE  NO.  OF  WORDS  TO  BE  SKIPPED  AT  THE  START  OF  TAPE10 
C 

READ(5,1000)  M, ITAP, NSKIP, MATZ 
N  =  M*(M+1) /2 
K1  =  N  +  1 
K2  =  N  +  M 
C 

C  READS  IN  ELEMENT  STIFFNESS  MATRIX  AS  A  LOWER  TRIANGLE  INTO  ARRAY  A 
C  READS  IN  SOLUTION  VECTOR  AND  TACKS  IT  ON  TO  THE  END  OF  ARRAY  A. 

C 

IF(ITAP.EQ.l)  GO  TO  50 
READ(5,1010)  (A(I) ,  1=1, N) 

READ(5,1010)  (A(I),I=K1,K2) 

GO  TO  60 

50  READ(IO)  (SKIP, 1=1, NSKIP) , (A(l) , 1=1, K2) 

60  CONTINUE 
C 

C  CREATES  SEPARATE  ARRAYS  AKE  AND  D  FOR  THE  ELEMENT  STIFFNESS  MATRIX 
C  AND  SOLUTION  VECTOR. 

C 

DO  100  1=1, M 
DO  100  J=1,T 
K  =  J  +  I*(I-l)/2 
AKE(I, J)  =  A(K) 

AKE(J.I)  =  AKE(I, J) 

100  CONTINUE 

DO  110  1=1, M 
J  =  M*(M+l)/2  +  I 
D(I)  =  A(J) 

110  CONTINUE 
C 

C  FINDS  THE  EIGENVALUES  AND  EIGENVECTORS  OF  STIFFNESS  MATRIX  AKE. 

C 

CALL  RS(50,M,AKE,EIG,MATZ,EIV,FV1,FV2,IERR) 

C  FINDS  SCALAR  PRODUCT  OF  EACH  EIGENVECTOR  WITH  ITSELF  TO  FIND  THE 
C  NORMALIZATION  FACTORS  AND  PLACES  THEM  IN  ARRAY  A. 

C 
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Table  9.1  ELMOD  Listing  (Continued) 


140 

130 


DO  130  J=1,M 
FACT  -  0.0 
DO  140  1=1,  M 

FACT  =  FACT  +  EIV(I, J)*E1V(I, J) 
CONTINUE 

A(J)  =  SQRT(FACT) 

CONTINUE 


SCALES  EIGENVECTOR  MATRIX  EIV  TO  MAKE  IT  ORTHONORMAL. 


DO  150  J=1,M 
DO  150  1=1, M 
EIV(I, J)  =  EIV(I,J)*A(J) 
150  CONTINUE 


THE  VECTOR  OF  MODAL  COEFFICIENTS  IS  NOW  FORMED  AND  PLACED  IN  ALP. 
CALL  MATRIX(ATRANB,M,M,1,EIV,50,D,50,ALP,50) 

DETERMINES  THE  VECTOR  OF  MODAL  STRAIN  ENERGIES. 


UE  =  0. 

DO  160  1=1, M 

U(I)  =  0.5*ALP(I)*ALP(I)*EIG(I) 
UE  =  UE  +  U(I) 

160  CONTINUE 


PRINTS  OUT  THE  RESULTS. 


DO  170  J=1,M 
WRITE (6, 1020) 
WRITE(6,1030) 
WRITE (6, 1040) 
WRITE (6, 1050) 
170  CONTINUE 

WRITE (6, 1060) 
STOP 


J 

(EIV(I,J),I=1,M) 

EIG(J) 

U(J) 


UE 


1000  FORMAT (413) 

1010  FORMAT(6E12.6) 

1020  FORMAT ( 1H1 , / / / , 10X , 11HEIGENVECT0R , IX , I 3 , / ) 

1030  FORMAT (10X,E12 .6) 

1040  FORMAT (//,10X,31HTHE  CORRESPONDING  EIGENVALUE  IS,1X,E12.6) 
1050  FORMAT (//,10X, 2 6HTHE  MODAL  STRAIN  ENERGY  IS,1X,E12.6) 

1060  FORMAT (//,10X,26HTHE  TOTAL  STRAIN  ENERGY  IS.1X.E12.6) 

END 
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Figure  9. 1  Main  Pre-processing  Program 


Figure  9.3  (See  Figure  9.1 ) 
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Figure  9. 1 3  Nodal  Coordinates  and  Integration  Points 
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Figure  9.16  Shell  Unit  Intersection  Overlay 


I'igurc  9.1  7  Shell  Unit  Loads  Overlay 
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Figure  9.18  Element  Unit  Loads  Overlay 
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Figure  9.20  User  Subroutine  Loader 
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Figure  9.22  Data  Transfer  from  STAGS1 
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Figure  9.25  Eigenvalue  Solver 
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Figure  9.27  Stiffness  Matrix  Formulation 
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Figure  9.28  Element  Stiffness  Data 


Figure  9.30  Overlay  For  1st  Variation  of  Strain  Energy 
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Figure  9.31  1st  Variation  of  Strain  Energy 
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Figure  9.34  General  Ordinary  Differential  Equation  Solver 
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Figure  9.35  Plasticity  Overlay  (Dynamici 


Figure  9.36  Overlay  For  Solution  Strategy  And  Output  Control 


Figure  9.38  Overlay  For  Stress  And  Strain  Computation  -  I D  Elements 
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Figure  9.39  Overlay  For  Stress  And  Strain  Computation  -  2D  Elements 
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